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Electric Polarizability of a Short Right Circular 
Conducting Cylinder* 


T. T. Taylor' 


(April 29 1960 


\ method similar to that employed by Smythe [1]? for calculating the capacitance of 
a freely charged short right cireular conducting evlinder is used hoc eee errs the electric polar- 


izability tensor in the principal axis system for such a evlinder. Calculations to an accuracy 
of —— ttely five significant figures are carried out for evlinders with radius to half-length 
ratios of! .1,2,and 4. The results are applicable to the design of artificial dielectrics. 


1. Introduction 


A solid conducting object in the form of a short, that is noninfinite in length, right circular 
evlinder makes a suitable element for use in the construction of an artificial dielectric. For 
such an object, therefore, it becomes important to know the electric polarizability tensor 
a, Which relates the induced dipole moment p, to the inducing field /, according to the equation 


pi=ayl;, (1) 


with the evlinder situated in free space. Once a,, is known, the effective dielectric constant of 
a spatial array of identical and identically oriented evlinders can be calculated according to 
the methods of Kaprielian [2]. 

A centered coordinate system in which the z-axis coincides with the axis of rotational svm- 
metry obviously constitutes a system of principal axes for a short right circular evlinder. 
The components @,, and a,, are clearly equal and are denoted a@,, where ¢ stands for “trans- 
verse’; similarly, @.. is denoted a,, where / stands for “longitudinal.”’ The evlinder cross section, 
with nomenclature, is illustrated in figure 1. This article describes a method which has been 
used successfully for calculating @,, and a@,, for eylinders with radius to half-length (a/b) ratios 
of 4. 4. 1, 2, and 4.) Rationalized M-.K.S. units are used throughout. 

The problem which presents itself is that of finding, in the external space, a Laplacian 
function which has no tangential gradient at the surface of the evlinder and which reduces to 
the potential of a uniform field at infinity. With such a potential, there is associated physically 
a surface charge density on the evlinder proportional to the local normal gradient and a field- 
free condition in the interior. It is the presence of this density which gives rise to the induced 
dipole moment. Since there is no known coordinate system in which a evlinder of finite length 
with closed plane ends forms a coordinate surface and in which Laplace’s equation is separable, 
the present problem must be regarded as intractable from the point of view of conventional 
methods and another method involving an arbitrarily good approximation to the surface charge 
density must be invoked. 

As a mathematical abstraction, let the evlinder be regarded merely as a geometrical con- 
struct to which the surface charge distribution just mentioned is firmly affixed. If that portion 
of the potential which corresponds to the uniform (applied) field is now subtracted, the remain- 
ing potential will be that due to the charge alone and will be found to correspond to a localized 
external field with dipole-like appearance and a perfectly uniform internal field which is equal 


* This article is based upon a portion of a thesis which has been accepted by the faculty of the California Institute of Technology in partial 
fulfillment of the requirements for the degree of Doctor of Philosophy Invited paper 
University of California, Riverside 
Figures in brackets indicate the literature references at the end of this paper 
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FIGURE 1. Cross section of the cylinder with nomenclature 


and opposite to the applied field. The problem may therefore be restated in terms of finding, 
on a finite evlindrical surface with closed plane ends, a charge distribution which, when acting 
alone, generates a uniform field in the interior. The technique employed in the present article 
makes use of this viewpoint and consists in setting up, on the surface, a charge density function 
which is completely determined by a finite number of parameters, then solving for the values of 
these parameters so as to obtain maximal uniformity of the charge field or, alternatively, 
maximal cancellation of the applied field, within. This method is similar to that employed by 
Smythe for calculating the capacitance of a freely charged eylinder [1]. The research reported 
here was performed under the direction of Prof. Smythe, to whom the author is indebted 
for many valuable discussions. The cooperation of the personnel of the California Institute 
of Technology Computing Center is also gratefully acknowledged. 

A second article, now in preparation, applies this method to the problem of determining 
the magnetic polarizability tensor under the assumption of negligible field penetration. — If 
both tensors are known and if the wavelength is long compared with element spacing, it becomes 
possible to calculate the wave propagation constant of artifical media whose elements are cyl- 


inders of the type considered here. 


2. Functions for Expressing the Charge Densities 


A system of generalized orthogonal polynomials and weighting functions has been con- 
structed for the purpose of expressing the charge densities on the side and ends of the evlinder. 
These polynomials involve an argument wv, which represents either 2/6 or p/a as required, and a 
weighting function (1—w?)’, The parameter v is given the value of minus one-third in order 
that the charge densities may be asymptotic to a constant times /“' as /-+0, where / is the 
distance from the edge of the eylinder. The product of the polynomial of order m and the 
weighting function is denoted by the symbol y,, and is defined as follows: 


1)”"2-*P'(m+y+- 6) 
(1—u*)*tF(- m, m+ v+¥yte Cyt Fu). (2) 


YnK ‘,o,v,U) - 
Yl: Piy+Or(m+1+ yp) 
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If the parameter y is given the value of one-half, this function becomes suitable for use on the 
side of the evlinder; if unity, for use on the end. The parameter ¢ controls the parity of the 
function and is zero for even parity, unity for odd parity. Another function, involving the 


polynomial alone, is given by 


(—1)"2'**(2m+94+74+ OT (m+ 74+-y7-4 


ms 5 Ys U) = -\Pr Fy (— mm + ety Sy+ ou’). 3 
Ym hy+ OrCm+1) vii YT ST! (3) 


The multipliers in (2) and (3) are designed to simplify both the orthogonality relation and 
certain integral transforms (see app. A) which play an important part in the analysis. The 
former is simply 


ae | —_ 
VW, U2? dub mn. (4) 


The polynomials themselves are adaptations of Gegenbauer and Jacobi polynomials; expressions 


for y,, and y,, in terms of Jacobi polynomials of argument (2uw°—1) are possible. 


3. The Longitudinal Problem 


In the longitudinal problem, a uniform electric field of magnitude EF is applied in the 
positive z-direction. The assumed charge density on the side depends upon the coefficients rp 
and r,,, 0<m< N,—1, as follows: 


: SA - {1 l 2 } 
0. ==)? ( ~) ‘ ly ht Do 1 ndn( » @ 3? h ), S < ). (5 
e ft 
0 2 b. 
The linear term, 
2 —- /1 2 a 
rn palo ¥ (5; 1,0.) ), (6) 


is called the “basic” term and is included in order that (5) may more easily approximate the 

true distribution especially when a/b is small. 
The potential generated in the interior by the charge density (5) acting alone is easily 

expressed as an integral transform: 
V.(9,.2) a f°T2 £ . ae : ” 
“ : | | R(z’) sin ke’dz’ | Koka) (kp) sin kz dk. (7) 
ck ¢ Jo Le Jo 

The Fourier sine transform of the charge density, enclosed in square brackets above, may be 


found with the aid of appendix A. The function /,(kp) sin kz, which is regular in a neighbor- 
hood of the origin, may be represented as a spherical harmonic expansion about the origin. 


Thus (7) becomes: 


V.(7,0) ab (2\'? J3,.(kb) X=! (—1)"J,2 a v(kb) 
’ Sy ‘ 3/2 ee . 2m +3/2+ # . ¢)2Ptl . 
aD 21 (-) |, E (ebyu2 baa!" (Eby vet |x (ka) (ke dk 
(—1)? (r\2?" ) 
(2p- The, P2p+1 (COS 6). (8) 


This expansion converges within that sphere which has its center at the origin and which 
passes through the nearest point occupied by charge, in other words, for all ra. 

Consider next the ends of the evlinder. The assumed charge density on the upper end 
depends upon the coefficients ¢, and ¢,,, 0<m<N,—1, as follows: 


Nem1 l p 
2 +> tn¥m ( 10-39") p<a. 
€ pal e)= m=0 


(9) 


0 p> 4a, 
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The constant term, 
rt p 
t» hb o(1,0,0, , 
tw = ) (10) 


is now the “basic” term. Again, the potential generated in the interior of the cylinder by both 
the charge density (9) and an equal and opposite density situated on the lower end, is easily 


expressed as an integral transform: 


Por Oe ee oe. as . 
a. | [ T(p’) Jy(kp’)p’dp i " STi(kp) sinh ke dk. (11) 
ch c J 0 a J ; 


This, in turn, may also be given as an expansion in spherical harmonies which converges for 


all r<b: 


V(7,8) 2a’ {* Ji(ka) , Neb (—1) "Joma ar(ka) ' , | 2p+l 
a > a =e VIF (ke)? dk Co) Parts (cos), 
ch p=0 C Jo ka m=0 (ka) , (2p +1)l\e vail 


(12) 


Evidently the sum of (8) and (12) constitutes an expansion of the interior potential due 
to the total charge distribution. Although this expansion converges only within the largest 
sphere which can be inscribed in the evlinder, it can be continued [3, p. 196] to any interior 
point of the eylinder. This expansion of the interior potential may be expressed in the follow- 


ing form: 


V(r,0) = "sb —_ =m red ~' mem p\rrr! > 
Au EB ry p + te lm Zz Pp + tx p -+ > & Ba NX p ( ) / op+1 (Cos A), (13 ) 
ck p=0 m=0 m=0 a 4 


where the .V’s are obtained from (8) and (12) by integrating [4, p. 137] over &. These ’s 
become the matrix elements in a system of simultaneous linear equations which may be solved 
for the rp, 7m, t, and ¢,, such that the coefficient of (r/c)??*"P.,,,(cos 6) will be unity for p=0 
and zero for as many p>O as possible. A solution of this type approaches exactness as the 
number of simultaneous equations approaches infinity; however, good accuracy was obtained 
with only eighteen equations. Among these is included an additional relationship known as 


the “edge condition,” 


N, =! a l ; N, cig 
Dra—(z) 35 tn=0, (14) 


which insures that the side and end densities match one another as they both tend to infinity 
at the edge. The expressions for the matrix elements are as follows: 


(—1)?+™ 22» T (om + p- 5 )E (om + p+ >) 


Er 


: 3 5 b? 
“aFi( me P +°, m—p+1+ 4; 2m-4 ~-+ v; x); (15) 


2 C 


x" | 5 
r ( \r (2p+2)P (2m+5+y) 


» 


| aoe 3\ 
(—1)™2°?-"r ( m+p- (m+ p+1) 
Pr35 / 


r(5) Pep+2) r(2m+2+¥) 


(Xe) 


, 3 . 
-aP( m + pts m—p-+- 1+; 2m+2--7; a) (16) 
_ Cc 


'¢ m 
«lp 


" . '" > ah reb ° ° 
Separate expressions for 7}, and 4, are not necessary since these elements may be obtained 


merely by setting »=0 in (15) and (16). 
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The dipole moment in the longitudinal problem has only a 2-component and is found by 
integrating zo over the entire surface of the cylinder. The orthogonality properties of the ¥,, 
and y,, greatly simplify this integration, which results in the following formula for the longi- 
tudinal polarizability : 


a, 2b 1"(3/2) | ; 
, : ” 
reo Ee aE 2°51 (13/6) }e[ oT a *1'(5/3) toh (17) 


where / is the geometrical volume of the cylinder. 
4. The Transverse Problem 


In the transverse problem, a uniform electric field of magnitude £ is applied in the positive 
y-direction. The assumed charge density on the side depends upon the coefficients s, and s,,, 
0<m<N,—1, as follows: 


8-1 S\ Yn 4 0; —> ;)p 2|<b. 


cos @| = (18) 





icos 


o. 
ey it 
The potential of this density in the interior of the cylinder is given by the integral trans- 


form 


\ ra “|. [2 [" S(2’) cos kz’ dz "| Kila) 1 (kp) cos ke dk cos $. (19) 
c 


Its spherical harmonic representation, convergent for r<ca, becomes 


V'.(r,8.0) > es €- 4 { [ J, 2(ko) S's, 1) Jen1/04.0(ke) 
ch p <6 ¢ Uk)! 8 " (kb)! rs 


. +1 nd) ol Lh 1 ‘ 
-K, (ka \(ke)*? dk (2p >) ( -) / 2p 4 1) (COS 6) cos d. 20) 


The assumed charge density on the upper end involves the coefficients w, and w,, 0<m 
<N,—1, and is given by: 


N.—1 
p — = l p 
ww, Su Vn l,l,—-s ); psa 
={cos @|W(p) =|cos ¢] o 6 o¢ (21) 
0, pa. 


The potential generated in the interior of the cylinder by both this charge density and an equal 


density on the lower end becomes 
VA , \ tae 
2 adhe “YI, fon “(pS (kp’)p’dp’ |e S, (kp) cosh kzdk cos ¢. (22) 


Again, the spherical harmonie expansion, convergent for rb, is given by 


V(r,6.0) a | a rs) (—1)"Joma240(ka) 
’ Wm ; 
ck + 24 (ka)! 


l ~~ : 
-e~ (ke)??* dk ( 91 ( ~) P34, (cos @) cos @. (23) 
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Once more, the sum of (20) and (23) constitutes an expansion of the interior potential 
due to the total charge distribution, which expansion converges within the largest sphere that 


can be inscribed in the cylinder. The expansion may be expressed as follows: 


V(7,6,0) 


N,-1 N.-1 \2p+l 
7 
™ "sh a) "sm . Yer > "em | _ a ‘ 
cK : > s,) > + S Sm >» TU ») at ST a = > ] ( - ) / op+1(COS #) cos ¢, (24) 
4 p=v0 m=0 m= 


where the Y’s are again obtained by integrating over k. It is desired that sy, 


Sm, Wy, and w,, 
— > mea l " ia 
be such that the coefficient of (r/¢e)?? T. (cos 6) cos @ be unity for p 


OQand zero for as many 
p>0 as possible, and a system of simultaneous linear equations must be solved. 


The edge 
condition now becomes 


N,-1 a 1/3 N,-1 
— =I _— 
DS &a—(z) 2) wa=0, (25) 
m= ) m= 


and the expressions for the matrix elements are: 


(—1)?+mQ2p T( m- p+3) r( ia p+s5) 
Y=. : 2 


; r(3) rep 30 (2m+5+») 


2 h\2m+2 
(Ge) (e) 


to 


’ » 9 | b m oe 
oF, ( m+ p+s m—p+1- y,2MTSTY; > ); (26) 


(—1)"2?-*P (m4 p+s)ri m+ p+2) ve . ; 
Snape (2)(2) (mint m—pitsnam sine) 


ce 
+ 3)P(2m+3-+ pv) 


Note that is and ty are obtained by setting vy=0 in the above. 


The dipole moment now has only an s component and is found by integrating so over 


the entire surface of the evlinder. The expression for the transverse polarizability is: 


a (1/2) al | l 
at —_ Sp T ans So + Ww), T , , ‘ We = (28) 
My€) 2?7/°T' (7/6) bi 4 2°/°T' (8/3) 


5. Results 


The two distinct tensor components of the electric polarizability are given in table 1. 


TABLE | 
alb ai Poe are Coe 
0 20 2 0000 
ly 15. 071 2.3151 
ly 7. 0066 2.6115 
| 3. 8614 3. 1707 
2 2. 4325 4.2173 
t 1. 7507 6.1814 
x L. ooo 

t=2rath 


These results may be compared with the polarizability of a 
to Sze In rationalized M.K.S. units. Actual values of 
cluded in appendix B. 


conducting sphere which is equal 
the charge coefficients are in- 
An estimate of the accuracy of the results was obtained by solving 
the problem repeatedly with larger and larger numbers of equations and observing the limit 
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toward which the calculated value of the polarizability tended. It was found that, for a maxi- 
mum number of equations equal to 18, the polarizability could be found to 5 significant figures 
with some doubt about the least significant figure. As might be expected, the accuracy is best 
for a/b equal to unity, in which case the least significant figure is in error by at most one unit. 

As an additional check upon the calculations, the potential, V, at the pole and the electric 
field, /., at the equator due to the charge alone in the longitudinal problem were calculated 
and compared with the corresponding values associated with the applied field. The small 
errors involved can, by appealing to the appropriate gradients, be translated into local dis- 
crepancies, Ap at the equator and Az at the pole, between the actual surface and an imaginary 
surface which is an equipotential under the joint influence of charge and applied field. A 
corresponding check, in which the potential, V, at the equator and the electric field, Z,, at the 
pole due to the charge alone were compared with the corresponding values associated with the 
applied field, was carried out for the transverse problem. The results of both checks are given 
in appendix B and are seen to be very satisfactory. 


6. Appendix A: Hankel Transforms of the /,, Functions 


It will be shown that all of the integral transforms of the ¥,, functions used in this article 
are special cases of the general transform: 
*! my ov 
| ; (—1)"P(m+1+-»)2°Fom414040(t) 
ur (1—u?)" PY? (2u—1) J, (tu)udu= \ = - = : o\ : 
P(m-+-1)t** 


Re o >—1; Rev >—1; t real and >0, (Al) 


0 


which will now be proved. P)?"' (2u?—1) is a Jacoby poiynomial [5, p. 168]. Denote the in- 
tegral in eq (Al) by f(t). Then 


. aa r ))\2r+e "= 
s(o)=>>5 (—)"(t/2) 1) | u??+2o(] — uy?) PP? (2u?—1) udu. (A2) 


fa (r+ (r+e4 


Letting 2u?—1=v, in [4, p. 284], 


ad | , , ‘ 
, r(/+1)0(m+14+~)0(/—o+1) 
a" 1—v? PP’ (202®—1)udu= 
" , 27 (m+ 1)P(l—e— m+ 1) (v+-/+ mf 


- ; Rel >—1;Revy>—1. 
2) 


(A3) 
Let /=r+o; Reo -]; 
: I'(m-+-1+-v) (—1)’(t/2)**¢ 
f(t)= : ; Reo >—1; Re »>—1. A4) 
y 20 (m+1) FU CTr— m+ lr(v+r+oe4+ m2) ( 
Make a change in the index of summation; let r=m-+-k: 
-1)"P(m+1+7) & (—1)*(t/2)*+2"+¢ : 
f(t) a ; rn - — . (A5) 
7 21 (m-+1) = FU+k)P(2m4+2+04+r+k) 
This is: 
. (--1)"P(m+1+4+ )/t\'-’ 
f(t)- (=) — Samsisero(t), AG 
mn 20(m+1) 2 matisoes\ — 


and the transform of eq (A1) is true. 
From [5, p. 170], one finds that: 


(—1)"%(m+1+ce) ,, ‘ ” 
Pe” (2u?—1) 'm+-1)P(o+1) of (— m,m+vt+at+1;¢+ 1; u*)- (A7) 


141 











If o is set equal to y+¢—1, the transform becomes 


‘ m 
on (—1) Jom rtar(t) 
Wnl¥ Cv Udys¢ (tu) (tu)rdu at ak; 


(AS) 
Now if the parameters y and ¢ are allowed to take all values relevant to the present problem, 


one obtains: 


Side, even: 


"] ] T are ¢ -1)* Jans rn4y(t) 
Yn ( 520,0,0) Cos tudu (7) ithe od (A9) 


Side, odd: 
"! l w\!/? (—1)" Joms3/240(t) 
| Ym (571m) sin tudu (5) ame ; (A10) 


ia] - a 
End, even: 


e) 
| Vn (1,O.v, a) J) (tuudu (All) 


End, odd: 
= 1 » ‘ 


api t) 
pits ; - Re »y>—1. (A12) 


*1 
Wn (lly, u)J,(tuudu 
J0 


Because of the special values taken by y and ¢, the singularity at the origin which might have 
occurred in (tu)%J,,¢-;(tu) does not appear. Both sides of eqs (A9) through (A12) become 
entire functions of ¢ and no restrictions need be placed upon f. 


7. Appendix B: Values of the Charge Coefficients 


The values of the charge coefficients and the results of the checking procedure at pole and 
equator are given in tables 2 through 6. The notation “.1(p)”? means “.1 times 10? ” 


TABLE 2 TABLE 3 
ajh='4 h=! 
— | 
| 
Longitudina! Transverse | Longitudinal Transverse 
m Tm Sm m ry 8» 
Basic +(), L2QO05187 (+0 +(), 40BS7H37 (+0 Basic 0. 19225150 j +0. SSYUH5487 (+0 
0 +0). 236032900 (+1 +1). L46007S87 (+1 0 +1). ISSB67T9IT (+1 +0). 16055212 (+1 
l —(), 113903441 (+0 0. 43230914 (+0 l 0). 488386001 l 0. 35500127 (+0 
2 —(), 19185823 (—1) 0. 65706908 (—1 2 0. 54930270 (—2 0. 438967119 (—1 
3 —(). 52174681 (—2 —(). 17221441 (—1 | 3 0). LOSS2670 (—2 0. LOLI9125 (—1 
4 —(). 1484165 (—2) —(). 54219231 (—2 i (), 28834773 (—3 0. 28016839 (—2 
5 —0), 53665621 (—3 0. 17807221 2 ) 0. 54147255 4 0. 7TS539136 (—3 
—). 16813125 (—3 —(), 5407302 (—3 6 0. 11308514 4 ), QIS71OL (—3 
7 —(), 48283874 (—4 —0. IHBASHOY (—3 7 0. 20056279 (—5 (). 43956298 (—4 
. —(). 12195426 (—4 0. 41591873 (—4 s —(). 27611488 (—¢ 0. 74296836 (—S5 
q —), 25981007 (—5 0. S8O90900 (—5 4 0. 25757920 (—7 0). S53901338 6 
10 0. 12042175 (—s 0. 49561337 (—7 
10 —(). 44409406 (—6 —(), 15241120 | 
11 —). 6727317 (—7) —(). 19470442 (—6 Equator E/E 1. GOOO0003 ' ak cos @= +1. 00000009 
12 —(). 47918712 (—S —(). 16423381 (—7 check 
13 —(), 20011235 (—¥ —(), 68402739 4 - 
Ap/a Negligible Negligible 
Equator E,/ E= —1. 00000002 'Vak cos é= +1. 0000002 
check | m tm Wm 
Ap/a Negligible Negligible Basic +(). 46038430 l +41), SUNBSTOS » 
m te Wm | 0 +4), 233570 (+1 +0). 13797193 (+1 
| ] 0. 82663114 ] +1). S840603 1 l 
Basic +4). 15200248 (—1 0. 17172000 (+0 2 (0). 76227915 (—2 +4), 204651 (—1 
| 3 0. 12692924 (—2 +1). 7724312 2 
) +(). 35624820 (+1 +). 14602027 (+1 { 0. 165849013 (—3 +), 15415744 (—2 
l —(), 39128623 (—1 +4), 32011553 l 
Pole check hk +4). WUU4UUUT Kp' KE cos @ 1. OOO02004 
Pole check hE = +0. Y9999875 Ep/E cos @ 0, 99902328 
Azih Negligible +1). QOOOO890) 
Az/b 0. (OO00037 0. OO022077 
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TABLE 





th=) 
Longitudinal 
m rT 
Basic 0. STOLOS1L7 l 
rT) +), 149210868 (+1 
1 0. 83747159 (—2 
2 +(), 230440N7 (—2 
; +() LB471808 (—2 
{ +() 48100777 (—3 
+1). 146008304 (—3 
t ‘ 4) 
7 +4). 34949152 (—5 
Equator EWJk 1. O@O00001 
check 
Ap/a Negligible 
m f 
Basic +0. 87113326 (—1 
0 +0). 16193448 (+1 
1 0. 11224086 (+0 
2 0. 14340444 (—1 
3 0. 34812949 (—2 
4 0. YSSTOTHL (—3 
5 0. 23486408 (—3) 
6 0). 44962281 (—4 
7 0. 45407822 (—5) 
Pole check hk = +1. 00000001 
Ac/h Negligible 
TABLE 
a/h=2 
Longitudinal 
m I 
Basic O.H5161515 (—2 
0 +4), 12482160 (+1 
1 +0). 24818534 (—1 
2 +0, VIOGOIST (—2 
3 +4), VSTOTYOL 2 
{ +0), SOT4ISOOO (—3 
Equator E../F 1. 00000271 
check 
Apia +4), QOOOOT 21 
m f, 
Basi HO. 1HORSTHY (+0 
0 +4). LIGS3321 (+1 
l 0. 12442481 (+4) 
2 O. 171506381 (—1 
3 0. 42141901 (—2 
1 0. 12OUSH44 ( ») 
5 0. 34737952 (—3 
th 0. WOO47300 (—4 
7 0. WL11927 (—4 
s 0. 34450642 (—5 
4 0. 40094338 (—6 
1 0), 234533528 (—7 
Pole check LT +1. OOOOO002 
Azih Negligible 
23 TO oO 2 


} 


5 


Transverse 


+), 25340314 (+0 


+0). LWSTIOIT ¢ 
(). 32219422 ( 
0. 38977600 (—1) 
—(). VS990090 (—2) 





{ 
{ 
—(). 12162489 ( 

0. 12244613 ( 


'Vak cos = +0. 99999999 


Negligible 


Wm 
0. 13425511 (+0) 
+1), L47U2105 (+1) 
+1) SS140061 (—1) 
+0. 17095007 (—1) 
+1). SYTT4ASYS (—2) 
+4), 2460043 (—2) 
+4), 50421158 (—3) 
+(), 12420137 (—3) 
+1). 13689068 (—4) 
Ep/ FE cos = 0. Gaga 


Negligible 


Transverse 


+ 


1611166 (+0) 
+4). 23744087 (+1) 

0. 28135556 (+0) 
—(). 33694020 (—1) 
72222000 (—2) 


0. 11581174 (—2) 


Viak cos @ 


+4), GUUO9ONS 


—(). QOOO0005 


0). 22721442 (+0 
+1) ISG58501 (+1 
+4), 426935906 (—1 
+0). 12814404 (—1 
+0). 43243348 (—2 
+0). LA99S6H39 (—2) 
+4). 49000507 (—3) 
+1). 14126809 (—3 


+(° 330607692 (—4 






+4). GI418057 (—5 
+1). 74203582 (—6) 
+4), 46182535 (—7 


Ep! FE cos @ 1. QOOO0003 


Negligible 








TABLE 6 


Longitudinal Transverse 








m Tm Sm 

Basic 0. 91407486 (—1) +0). 36081774 (+0) 
0 +. 11655470 (+1) +. 27748758 (+1) 
l +. 92149445 (—2) 13470990 (+0 

Equator ko k= —0, 99979406 lak cos é=+1. 00001613 

check 

Ap/a 0. OOOO 5287 +). 00000530 

m tan Um 

Basic +(), 1GS31628 (+4) —(). 10545599 (+1) 
0 +1). 92028755 (+0) +4). 15388757 (+1) 
1 0 +4). 82057625 (—1) 
2 —(), 2 +1). 25448502 (—1) 
3 —() +). QOS57682 (—2) 
4 0). 23904452 (—2) +(). 30007948 (—2) 
5 0. 81543454 (—3) +4). 15309282 (—2) 
ti) —(), 26584526 (—3) +1}, 55635828 (—3) 
7 —{(), TRRYBO56H (—4) +4). 18020312 (—3) 
s 0. 2466149 (—4) +4). KI2H2T7O8 (—4) 
if) (0). 44553523 (—5) +1). L1G35018 (—4) 
10 0. 77488570 (—6) +(), 21335661 (—5) 
1! 0. 10037598 (—6) +(). 28045170 (—f6) 
12 0. 85747540 (—8) +4). 25757752 (—7) 
13 0. 36132625 (—9) +). 11256670 (—8) 

Pole check lhe =+1. 00000001 Ep/ F cos @= —1. 00000001 
Az/t Negligible Negligible 
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Distribution of Quantiles in Samples from a 
Bivariate Population* 


M. M. 


(March 3, 


Let F(x,y) be the joint ay yey function of 


and Fa(y) the 


a quantile of F(x) 

er is drawn and t 
Vi Nj and Y,< Y; if i<j. Let i and 7 bet 
jin <F.(B), and let M be the number of eleme 
joint 


function f r,y). Let F(x) be 


spectively. Let a@ be 


(X,, Y,), A 1,2 


n, 


distribution of (AJ,N;,} is obtained 
estimates and confidence limits on the parar 


1. Summary 


The exact distribution of quantiles as well as its 
limiting form is well known in the univariate 
Mood ? investigated the joint distribution of madiae 
in samples from a multivariate population and 
showed that it is asymptotically normal. In this 
paper the exact joint distribution of qui antiles and an 
auxiliary statistic in samples from a bivariate popu- 
lation are obtained and it is shown that the joint 
distribution asymptotically trivariate normal. 
The auxiliary statistic utilized to estimate the 
correlation coefficient between quantiles and also for 
setting up confidence limits on it. Estimates for an 
ordinate of a univariate probability density function 
pdf) are also obtained. These are used in setting 
confidence limits on the quantiles. 


“ase 


Is 


Is 


2. Introduction 


Let F(v,y) be the distribution function of CY,Y) 
a pdf f(r,y). Let the marginal distribution 
function of Y be denoted by F\(2) and that of Y by 
Fy). Let @ be a quantile of Fy(7) and B be a 
quantile of Fy(y). It is assumed that the first and 
second partial derivatives of F(r,y) are continuous 
ina neighborhood of (a8) and that f(a@,8)#0. A 
random sample CY,,¥4), A=1, 2... ., ”, is drawn 
from F'(r,y). Let the sample values of .Y be ordered 


possessing 


so that 

Xi NX; ie 
Similarly, let Vy Y emis eae All such samples 
for which NX'=N’, or V'=Y" for different indices 7 


and j may be exeluded from ‘consideration, as these | 
Let ¢ and 7 | 


samples form a set of probability zero. 
be the integers such that 


Un<Fy(a)< ((+-1)/n3 j/n SFB) G+ 1)/n. 
*Contribution from the 

Boulde r, Colorado. 
Harald Cramér, Mathematical methods of statistics, p. 368 (Princeton Univ. 

Press, Princeton, N.J., 1946 

*A. M. Mood, On the joint distribution of medians in samples 


‘variate population, Ann. Math. Stat. 12, 268 (1941), 


from a multi- 


and 6 be 


National Bureau of Standards Boulder Laboratories, | 


Siddiqui 
1960) 


(VY, ¥), possessing a probability density 
marginal distribution functions of X and Y re- 
a quantile of Fo(y). A random sample 
he values on each variate are ordered so that 
he greatest integers such that 7 n SF ila) and 
nts (VN, Y) such that Y< VN; and Y¥<_ Y.. The 
and is shown to be asymptotically normal. 
neters of interest are also given. 


Hence 7/n = F\(a) —6,/n, j/n=F,(B) —6,/n, 0 <8, 6g 1. 
If \f denotes the number of elements CY,Y) in the 
sample such that YY) and Y<Y%, then M is a 
| diserete-valued random variable with possible values 
0,1, , min (@—1, 7—1). 

First the exact distribution of (MZ, X4, Y%) is ob- 
tained for a fixed n. Then the asymptotic « distribu- 
tion of the standardized variates found when 

i,j so that 7/n—>F\(a@) and 7/n—>F,(B). 


Is 


>= 


“, 


3. Joint Distribution of the Three Variables 


Let us take a Euclidean plane (2, y) to represent 
the sample points. Given two numbers zy and yp, 
the lines r=s) and y=yp divide the plane into four 
quadrants, namely 


(); (x, y): tao, YC Yo}, 
()» (x, y): £20, Y>Yo}, 
()s (x, y): © >To, YN Yo 
td; (7, y): >to, Y >Yo 
Let p Pr{(X, Y)eQ,;}, t=1, 2, 3,4. Then 
pi=F (xo, Yo), Po= F(a) — F(x, Yo) 
Ps= F (yo) — F (xo, Yo); ps=1—F (20) 
— F (yo) + F(x, Yo). 
NX) will fall in (2, 2o+da) and ¥4 will fall in (y, 
| yot-dyo) in all samples CX,, ¥), #=1,. . . , n, such 
that the following statements hold simultaneously: 
| (1) ¢—1 points fall on the left of the line r=2p, 
n—? points on the right of it, and one point on it; 


(2) j—1 points fall below the line y n- 


points above it, and one point on it, 


Yo, J 


CX’ yo) 
Yo respectively. 


denote the points which 
rT . 
rhere are five 


Let (7,¥"’) and 
| fall on ray and y 
| distinct possibilities: 
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Case (1): X’<ap and ¥’< yo, 

Case (2): X’<sx) and Y’ >y%, 

Case (3): X’ >a) and Y’< yw, 

Case (4): X’>2x) and Y’ >y, 

Case (5): Y’= zo, in which case there is only 
one point (2 ,%), common to both 
the lines. 

We will consider Pr(M=m, 2 <X) <2 +dr, Yo 
< VF < yp tdyy) =P (m2, Yyo)drodyo. Uf Mm, then, 


in case (1), there will be *—m—2, j,—m—2 and 


n—i—j+m-+2 points in @, Ys, and Qs respectively. 
The probability of such a sample is given by 


P, (m,20,Yo) dred Yo 


n! m,,i-m~-2,,i ~m-2,,"n+m-i-i+2 of oF, iD 
Uptps “1! *» ite droll 
NU Or OY y 
m! (1— m—2)!(j7— m—2)! (n+ m—i1—j+2)! 
(3.1) 


where the partial derivatives are evaluated at (29,%). 
The contribution to P(m,r9,y) from other cases 
may be found similarly. The possible values of m 


in each case are determined by the rule that negative | 


factorials are not allowed. For example, in case (1), 


the maximum possible value of m is min(¢—2, j7—2) 
and the minimum possible value is max(0,7+j—1" 


—2). It will be assumed that the suffixes 1 and 2 
are so chosen that F\(@)<F(8) and hence i<j. 
Thus, with the convention that the terms involving 


where all the partial derivatives are evaluated at 
(19, Yo)- 


Hence the joint pdf of CY),¥%5), 


PA Lo, Yo), is given by 


i—1 
P(Lo.Yo) os P(m,25,Yo), Sto, Yom ™, (3.4) 
m ou 
and the probability distribution of AZ is given by 
P(m) | P (m2 ,.Yo) drydyo, 
m=0,1,...,%—1. (3.5) 


4. Asymptotic Distribution 


In this section the following will be found: 
(7) The limit of the pdf of (U,V), pwr) (the no- 
tation p is used generically to denote a pdf) where 


n/?(X; 
[Fi(a@) {1 


8) f.(B) . 
FB) }]' 


a) f(@) : a'/2()’, 
F(a) }|'? [F’(8) 51 


(77) The asymptotic probability distribution of 
= M/n as n>. 

The following operations are performed. 

(1) Write 


negative factorials be replaced with zeros, F(a,8)=F, Fi(@)=F,, FAB)=F, fil@)=hi, 
P (m,2,Y%o) =>) P(M,x),Yo), f, , i o*k 
: : (8B) =f, hh, ] 3 =], 
k= Or OY Or ov" 
m=0,1,...,i—-1, —©<m, yoo, (3.2) , 
‘ , ‘ : A / a /fo=Co, I,(a@) , Delf) Me 
where ?, is as given in (3.1) and Ordy * WNv=Cry  G2lJ2=C2, Ilk@)=fiy Ia'P)=Ss 
alae tes illite, alana th ( f, (tp) oF ) oF ; 
™ 2 } fiaj)= 
; 1 Pe Ps Ps Ji or ) dy 
P; 73 — — , 
m!(1—m—2)'(7—m—1)'(n+m—i—j+1)! 
ni pT p, "'pi-*-*pit* i+! =< f, (Yo) oF ) 
ul 2 { (Ue) — 
p _ va oz \"* ov 
, m!(1—m—1)!(j—m—2)!(n4+ m—i—j+1)! . (3.3) 
OF \/. oF 
yeaa — aiid 7 alin '( f, (2%)— f.(Y)— ) 
P pips") p’ fi a Mey dy 
: m!(i—m—1)!(j7—m—1)'(n4+ m—1—)7)! 
p. n! pt'ps m ‘ph m tp’ +m )+ f | Lo,Yo) 
*  m!(7—m—1)!(3—m—1)!(n+ m—1—j+1)P 7 
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whe 
(a,8 
of t 
[tis 
Als 


and 


Th 
ran 
ble 

()}- 


for 


wh 


Gl 


Hi 


ad at where all the partial derivatives are evaluated at 
(a8). It is observed that by our assumptions none 
of the quantities /,, fo, gi, ge, and f is equal to zero, | 
Itisfurther observed that F< F) since F'(a,8) < Fa, ). 


i by 
; Also, | 
3-4) TT flay) dy< Slay) dy=fi, 
and g2<fo, So that O<e¢,, <1. 
(2) For S>2, use Stirling’s approximation, 
S! y2a Sots 11+0(8 )}. 
(3.5) a rs 
(3) Write 
mon q, ] 4 dy. 
This operation is equivalent to replacing the discrete | 
| | re] | 
no- random variable M Ti by an continuous random Varla- | 
e ble Q@ over the range (0O,F,;—1/n). However, 
—M ni<o1/n. 
3) If by these operations P? (m,zo,yo)drodyy is trans- | 
I 9 formed into PAG Lo Yoddqd rod yo, then 
PAG. 2oYo) =(@( Qty Yo) "1 (q.ty.yo[14+Om-")|, (4.1) 
of 
| where 
GG, 2.Yo) 
a pips’ “py (py "+4 . 
q'( Fi —q)*'~-*(F,— 9g)? “(A — F, — F,4+- q)'- Pet 
(4.2) | 
| 
} Hi rye) he “pik q) “(FP qg) : oF oF | 
f' ] yf (Dar)? *(pops)?q' 2(] EF F, j q) 2 or oy 
, Of \ oF 
} -(1—F f+ aq)(F; "8 £,(2,)— ) 
pp K—- Ft QE @ (file) —9, oy 


: , Of. oF 
l pi . \( Ff l o- | 
+ poppy '\A—F,— F.+- q) (Fi —q) , ( fo(yo) > ) | 


where all the partial derivatives are evaluated at 
(YoYo). The terms in /7 arise severally from P;, Po, 
P,, and P, after the factor G" is taken out and are 
given in that order. The contribution from P. is of 
order 7 us compared to the other terms and hence 
absorbed in O(n>') in (4.1). 

(4) Make the transformations 


n'?(r,—a)f, n'* (yy —B) fo 


(FG—F)!? “TRF 


Ti 


so that 


\ F | l KF, )( l ~ ia 


dryly nf f ») dude 
J 1/2 


and p(qto,Yyo)drodydq goes into p(q,u,v)dqdude, 


| where p is used generically to denote a pdf and is not 


the same function from equation to equation. 
(5) Expand each function of (a,7) in Maclaurin’s 
series up to the terms of order n~', e.g., 


’ (Ip a )? Yo )* 
I= F+ (4o—@) 91+ (Yo—B) Gat = a+ =f 4 


+ (49—a@) (Yo—B)f + O(n-*”), 


where sy and y are replaced in terms of uw and v 
Then 


log G W(q) + (to—a@) (q Py kyf koe, — F(F,—F)(1 


FP) }- (Yo B)(q—F )ky fo kc. — F(F, 
F)(1—F)) win 


») 


(1 a)*k, f? kc? + (1 


— l 
2¢,) F(F,— F) (1 — F2) + (q—F)Ai} —5 (Yo 


B)*k, f3! kee + (1 —2¢,) F( PF, — F)A—F,) 
| (qY Fy) A, (Io a) (Yo -B)ki fifo hoe) Co 
FUR, —F)A—F\)—eF(F.— F) (1 —F,) 
+ FUR, — F)(F,— F) + (q—F)A;3} + O(n-*”). 
Here, ¢, and ¢, are defined in the beginning of this 
section ; A), ly, and Ay are some constants depending 
only on the parameters of F(z,y) and need not be 
specified, as it will be seen presently that (q—F) is 
of order n7~'/?; 
i/k, =F (F,—F)(F,—F)A\—F,—F,+ F), 
ko= FFA — F,— F,4-2F) — F?; 


and 


w)=4 lo : + (Fi —gq) log Path q) log pe 


I—-F.—F,4F 
i—y,—F+< 


» 


+ (l—F,\—F,4+-q) log 


It is easy to verify that w(q) has an absolute maxi- 
mum at g=F in the range of values of g and w(F) 

w'(F)=0. Expanding w(q) in a Taylor’s series 
nbout q F’ 


| , . 
w(q) = —shiko(g— F)*+Ol(q—F)'. 


147 





Each funetion of g is expanded about g=F and THEOREM. The asymptotic joint distribution of \ 
the transformation | the variates (.X{,Y;) is normal with parameters 
t= yk,kon(q—F) (4.3) EX;i=a, EY ;=8, U= 
is made with a similar transformation 7 on the Var F,(—F,) Vary” F,(1—F,) 
random variable Q. The pdf of (7,U,V), p(t.ur) — aa * 7 wa 
is then given by = 
‘oR it 
' ies I I F, wil 
9 3/9 » 1/2 l > ( oV Cia ,) > “s UW, 
p(t,u,r) = (29) ~-*/?(1—p?)~'? exp] —5 (t—uhy—vh,)’ : nf fe ali 
in 
u?-+2—2pur Integrating p(f,u,r) with respect to w and ¢, it js tal 
ee eee it +On-'*)) (4.4) | found that the distribution of 7 is asymptotically 
2 p | normal with mean zero and variance 
where | pa 
pO sie ff o}—1+h3+h34+ Qhyhop. 
VA, FL0—F)) 0 — W) 
Hence, from (4.5), @ Is asymptotically normally 
j vii lkoe, — F( FP, — F)(A— F,)| distributed with mean F and variance 
“= = —— 
y (1—p*) FL(1—F,) le 
oe (hy kon) 1-4 h? 1 hz T Zhihop). (4.8) as, 
{ 
Vk lkoco— F(F,—F)A—F;)| a on : , les 
h,= a = = = Me (4.5) The variate M/n, which is) discrete-valued, is ( 
v(1—p*) Fi (1—F) distinguished from @Q, which possesses a pdf.  How- 
ever, if q minh, 
The limits of variation of wand ¢ are (—o, @) and 
= a . s 
those of ¢ are (—d, yi, t dyn) where d, and dz | Pr (M/n<q) Pi Sq. (4.9) (lis 
are positive constants independent of nv. In fact, | 
a — 5. Confidence Limits on Quantiles 
d,=yk\k, max (FF, + F,—1), d,=~yky ko (F,—F). 
When n--o the joint distribution of (7/7,V) Since .Y) is asymptotically normal with mean a 
is a trivariate normal with mean vector (0,0,0) | and variance 
and covariance matrix V obtained from 
Zs ” - is 
of =F, (A—F;)/(nf}), 
l —h, -h, 
—— 21/1 22)-1 _ay-1 |. confidence limits ona with confidence coefficient 
' hb W+(—e") hyhz—o(1—p*) (1—y) are given by Nj+Z,o,;. Here, 0“\y<1, and 
—hy hyho—p(1—p*)~! h3+(1—p?)~! a 
(4.6) | _ thing ORS 
JZ, \2r 2 
Integrating out ¢, the asymptotic pdf of (77\V) 
is obtained as However, ¢, remains undetermined unless f, is known. 
Exact nonparametric confidence limits on quantiles 
‘) — (Qe) -1(1—p?) -"/? ex I :' are readily available.’ However, an alternate pro- 
p (u,r)= (24) ~" (1—p") weed tat Yh p’) . cedure is outlined here which provides an asymptotic Ni 
| estimate of f; and hence of o;. . ; id in 
Ly2—2our} | —o<u, r<o, (4.7)| < onsider the joint distribution of —(Xj+,Xs 9 
oni ” ; Mesa where as no, ?—o In such a way that C>. 
i/in—>F, while k/n and //n both tend to zero. k and 
so that p is the correlation coefficient between U7 and | / may be taken to be of order n*, O<a~ 1. The 
V and hence between X/ and Y/. The following | probability that t<.NXj<t+dt, t.< Ni. <t+dh, 
theorem may then be stated: tr<.N.,<t+dt, is given by p(tt),.t.)dtdtdt,, where 
é _ 2 T 
(tty. ts) n![ Fy (t,) |)" [Fy (0) — F(t) FOF (tz) — F(t) | F(t) |" A) A(t) hi (te) (5.1) OX 
) , 9©2 : . oo. 
' ahi (i—k—1) (k—1)'7—1) (n—7—-1)! or 


38.8, Wilks, Order statistics, Am. Math. Soe Bull. 44, 6 (1948); G,. E. Noether: 
On confidence limits for quantiles, Ann. Math. Stat. 19, 416 (1948 
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Mi of 


it is 


rally 
ally 


4.8) 


OwW- 


4.9) 


ant 
nd 


he 
ty, 
re 


Make the following transformations: 


__yn(X; a) fy U, Tt] 


= CY; 
y Fi 1—F) k 


ae 


with corresponding transformations from (¢,f;,f2) to 
Use Stirling’s approximations for factor- 
ials involving 7 and. Expand each function of u 
in Maclaurin’s series and let n->e. We finally ob- 
tain the asymptotic joint pdf of (0'07,,07,) given by 


UU, ) 


‘ u . (kf,)’ 1 kf (lf,)' l f 
) Us) u} 7141 —— Us “2 
= 2x (k—1)! iy * 
r<ucw, O<u, Wom. (5.3) 


The following statements are true asymptotically. 

(1) (07,072) form an independent set of variates. 

(2) 2kfl",=2nf, (Vi —Xj_,) is a X%* variate with 2k 
degrees of freedom. 

(3) 2/ ft’. 2Qnf (Nis. 
degrees of freedom. 


(4) Henee, 


Ni) is a X* variate with 2/ 


2n (X54. —Ni_wdhi 


variate with 2(k+4/) degrees of freedom, 


rr 
i- 


is a X" 
distributed independently of 
(5) Thus, 


, U(k+1) (V'—a)(k+/) 
Ww. . : ~ - 
f, (kl itll 2) yuk SL? A ) 
(5.4) 
is distributed as the ratio of a N(O0,1) variate and 





» 


p=—4F-1 kp '=F?(1/2—F)? 
2¢,—1 2c.—1 
h, o 9 hy ee =9 
\ l \ | - p- 
Fi—2F) l : : 
2 2 ot 2 
7 aT) T1én (20 


Now 0<¢,, @<1, O<F<1/2, 


C.—1)*-+ 


an independent xX? variate divided by its number of 
degrees of freedom 2(k+-/). 

W is independent of the parameter f,; however, 
the distribution function of such a ratio has not been 
tabulated, hence, at present, it is of little practical 


use. However, 1/f,; can be estimated by 


U,=(kU, 4+) /(k+D, (5.5) 
which has mean value equal to 1/f; and variance 
(k+/)~'fy?. Thus 
Sj V Fy ] — Kyl I. \ n (5.6) 
is an unbiased estimator of ¢; with known distribu- 
tion, still in the asymptotic sense. 
In practice, k and / may be taken of order n'’”. 


6. Confidence Limits on F and p 


The asymptotic maximum likelihood estimator of 
F is M/n which has asymptotic variance o3, given 
by (4.8). Solution of the inequalities 


—L4,09<Q—F<Z,0¢ 
gives 


h(Q) <F<h(Q), 


to obtain asymptotic confidence limits on F' with 
confidence coefficient (1—y). For large n, og may 
be replaced by Sg, where SZ is obtained from the 
right-hand side of (4.8) by replacing F by Q. Thus, 
to order n7'?, 100(1—y) percent confidence limits 
for F are G+ Z,So. 

Here, again, there is the difficulty of nuisance 
parameters ¢, and ¢. To show how to overcome 
this difficulty, the case of medians is presented. 


When F, = F,= 1/2, 


, (6.1) 





2p(2¢,—1) (2¢e,—1)}. 


and the expression | confidence limits on F' with confidence coefficient at 


in the square brackets has a maximum equal to | least equal to 1—y. 


2+2) =max (SF, 4(1—2F)} for variations of ¢, and 
® Thus, 
F(1—2F) , max(2F, 1—2F) 
a2, < 4 , <1/(4n). (6.2) 
de 2n r 4n _ ( 
The last bound is the least upper bound of the 


expression in the middle and is attained when F'=0 
or 1/2. 


Henee V+ Z,/(2yn) are asymptotic nonparametric 





Going back to the general case, it is observed that 


Q—F,F; 


(6.3) 
\ BFA —Fy)Q — F) 


is the asymptotic maximum likelihood estimator of p. 
Confidence limits on p are easily set up as p 1s a 
linear function of F, and r of Q. 
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A special case. If F(x,y) isa bivariate normal with 
correlation parameter p* then, for j= F,=12, 


C) Co ] 2 h, hy (), 
l l , ee 
I a""_ cos‘? O<cos”'p*<z, 
p*=sin5p=sin, (4/—1), 
» F(i—2F) 
O07 ‘ 
” 2n 


and 100(1—y) percent asymptotic confidence limits 
for p* are 


sin 3 yv8U(1—2@) /n } 


as Sin @ is a monotonically increasing function of # in 


(—w/2, w/2). 


4—14 


7. Generalization 


The results of this paper can be generalized to the 
case of a set of quantiles in samples from a multi- 
variate population possessing a pdf and satisfying 
certain continuity conditions. From the discussion 
of Mood (see footnote 2) of the joint distribution of 
medians in samples from a multivariate population, 
and from considerations of the generalizations of 
eq (3.3), it is seen that the asymptotic distribution 
of quantiles will be multivariate normal. To be 
more explicit, the joint distribution of a set of quan- 
tiles is obtained as a weighted sum of multinomials 
and multinomials are known to tend to normal under 
suitable linear transformations on the variates. Thus 
to specify the asymptotic distribution only asymptotic 
means and covariances of the variates are needed. 
However, these can be determined by considering 
the marginal bivariate distributions. 


Then the following conjecture may be stated: 


Consecture. Let F(x, . . ., 2,) be the distriby. 
tion funetion of CY,, ayy Boece ssing a pdf 
f(a, .» Zp). Let the marginal distribution of X 
be denoted by Fi(r) and that of CY,¥,) by F’s,(2,y) 


and corresponding pdf's by fi(v) and fy(2,y). Let 
Guns Ge . « «Xt, d ie Soar Pp 
be r; quantiles of Fy(r). Continuity of the 


derivatives of F up order p aut the 
(ay), ) is assumed for all gigs 


partial 
points 
Values 


to 
A» 
of ()), . - »jp)- Leta sample of size n be drawy 
and the vi alues of each variate be ordered so that 
N’(11)<.N’(22)- NX’(in), *=1,2,...p. 
Let Bua, } me! 3 o fe 4 ‘2 . e he 
the integers such th: at 
k,, n< Fila y< (k l)/n. 
the 
X 


‘ai 
Then 


Varinites 


jomet distribution of the 


is normal with parameters 


asymptotic 
"(ik;;) 
BX’ ( 


tk .,;) = ai; gun}, 2, Ts i=1,2,...,p; 


I (Qi; ;,Qlip,) 


nt (a 


ia (a F(a, 


)f (a@,,,,) 


Cov {X’(ak,,), 


Here Fy, (aj), im 


Ing Way: 


The author acknowledges the helpful criticism and 
comments by Dr. Edwin L. Crow on an earlier draft 
of this paper. 
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URNAL OF RESEARCH of the National Bureau of Standards—B. Mathematics and Mathematical Physics 
Vol. 64B, No. 3, July-September 1960 


Split Runge-Kutta Method for Simultaneous 


Equations 
John R. Rice’ 


(March 31, 1960 


Consider two simultaneous first order differential equations 2'(t)=F(a,y,0, y’ (bd) 
G(r,y,f). Runge-Kutta type integration methods are developed which allow different inte- 
gration steps to be used for these equations. These methods retain the desirable properties 
of Runge-Kutta methods, namely the self-starting property and ease of change of integration 
step. Two different approaches are considered and extensive experimental work is reported 
upon. Experiments are done both in situations where these methods are advantageous 
and where they are not. It is seen that these methods are more efficient than the normal 
{unge-Kutta methods if they are at all applicable and in ideal situations they give the same 
accuracy With 90 percent less computation. These methods are applicable to six degree 
of freedom missile simulations, for which they were developed. 


1. Introduction 


Consider two simultaneous first order differential equations: 


a ki t) ty) (1.1) 
P Zz .3 . A i) Ig . 

dt ry r ( J 

dA) / 

7 G(r,y,t) y(to) = Yo (1.2) 


where y(t) does not depend strongly on x(f) or varies much more rapidly than 2(f). In a 
normal numerical integration method for these equations, the integration step h must be chosen 
small enough to adequately integrate both (1.1) and (1.2). In this paper Runge-Kutta type 
methods are described which allow different integration steps to be used for these equations, 
These methods retain the desirable properties of Runge-Kutta methods, namely the  self- 
starting property and the ability to change the integration step easily.? 

The problem is defined in detail and two different approches to the development of the 
formulas are given in section 2. The analysis for third order integration formulas is given in 
sections 3, 4,5, and 6. In section 7 the results are stated without derivation for fourth order 
integration, 

Three systems of differential equations have been solved using these formulas with varying 
values of the parameters. The first equation is of the type suited for these formulas and they 
result in a considerable saving in computation. The results are discussed in detail in section 9. 
They point out that the second approach gives formulas which are considerably more accurate 
than the first approach— although one would not expect this beforehand. The second equation 
is of a type not suited for these formulas. The results are discussed in section 10. The third 
equation is of the type suited for these formulas except that a very high frequency low-amplitude 
oscillation has been added to x(f). The experimental results are somewhat erratic. 

In the final section a detailed discussion is given for the situations where these formulas 
are most useful and also comparison is made of their efficiency with that of the usual methods. 
It is seen that these formulas are more efficient than the normal Runge-Kutta methods if they 
are at all applicable and that in ideal situations they may give the same accuracy as normal 
Runge-Kutta with 90 percent less computation. One area of application is to six degree of 
freedom missile simulations, for which these formulas were originally derived. 


Chis work was done while the author was at Autoneties, Inc. and at the National Bureau of Standards as an NRC-N BS Research Associate 
For an account of the basic properties of Runge-Kutta see F. B. Hildebrand, Introduction to Numerical Analysis (McGraw-Hill Book Co., 


Inc., New York, N.Y., 1956 
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2. Problem Definition 


h denotes the integration step for (1.2), and Ah, where A is an integer, is the integration 
step for (1.1). Let ¢, denote t)+nh and let x, and y, denote the numerical solutions of (1.1 
and (1.2), respectively at t=f,.Round-off error is not considered in this paper. 

Consider two f-axes, 


_ Kh 
t, K t r IA t A 
- h 7 
Cua tmK-+1 eos t m+DK l f m+1)K+1 eS / m+2)K 


the first for (1.1) and the second for (1.2). It is now desired to obtain Runge-Kutta type 


integration formulas that integrate (1.1) in steps of AA and (1.2) in steps of h. 
(1.1) can be integrated from fxg tO finco« by a normal Runge-Kutta method. For a 
third order method the pertinent equations are: 


ky =KhF (12, Yn, ty) 7 
ki =KhF (2,+aWho, Yn tyiho, ta tyKh) 

ky = KhF (2,4 yak + (2-3) ko, Yn t+ Yahi + (ve —-Y3) ho, tr t+ yokth) 
hyo=KhG (4, Yn, ty) > (2.1) 
h,=KhG (2, +y1ho. Yn t+ yiho, th ty Kh) 


h. KAG (sr, t v3k) t (Yeo —¥3)ko, Yn 3h, T\Y2 3 )ho, t,t yoIvh) 





DP ¢m+1)K DimK T Boko T 3,h, T Bok». J 


The integration parameters ¥;, Y2, Ys, 8, 6, and 8, may be those of any third order Runge- 
Kutta method. Note that A, need not be computed for this integration. 

The main difficulty in integrating (1.2) is to obtain values of r(f) at the integration points 
between tnx and toasnx- A natural way to obtain these values is to extrapolate «(f) from f,,. 
The Runge-Kutta method is itself an extrapolation process and one extrapolation has been 
made in the integration of (1.1) from ty tO fonsiyyx. The values of ko, ky, and sy from (2.1 
may also be used to extrapolate r(f) from ¢,,, to the intermediate points. Let v,,<,, denote the 
extrapolated value of x(f) at free, ISySAK-1. Then, for the appropriate coefficients \,(/), 


3 
Ln+7 > > ADA 1. 
i=1 


Other estimates of r(f) are needed and they are obtained in the same manner. 
For a third order method the equations for integrating (1.2) from ¢,,¢.; tO fags.) are: 


dy(j) ha (mk + Ymk + rt mk + j) 


dy (G)=RG (ante Ds (G) iat s+ Maloy md uh ) 
i=4 


d,(j) NG ( tne + DD A (pk; 72Y mK +57 us, + (we ~ 3) dot mx + Loh ) 





YmK+j+1 YmK+j77 aly (7) T ad, (7) T aryl, (7). ! 
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Two basic approaches are considered for obtaining the integration parameters 4), po, Ms, 
a, a, and a, and the extrapolation parameters A,(j), . . . , Ay). The first is to consider the 
truncation error of the resulting integration formula for y(t). The parameters can be chosen 
so that the method is third order. This is done in sections 3, 4, and 5. The second approach 
is to determine the extrapolation parameters so as to make the extrapolation as accurate as 
possible and then to determine the integration parameters independently. This is done in 


section 6. 
3. Preliminary Computations 


For simplicity 27 2°\,(j)ky_, will be denoted by 3,. Likewise, 


A,(j): An(J) + Ans) +Anse(J) 
r,(j): ViAnsi (i) +YaAnse(J) 
A,(j): V7 Angi) +Y72Ang2(J) 
dn(j): V17V3An42(J) 


These expressions will occur quite often. The argument 7 will be omitted if it leads to no 
confusion. 

In order to consider the truncation error of the integration formula, various terms will be 
expanded in Taylor’s series. There will be some derivatives evaluated at (2nx«,Ymx,tm«) and 
some at) (ngs) Ynksjtmes)) Expansions will be obtained in this section that will allow a 
comparison of such terms. A Taylor’s series expansion of Z, will also be obtained. 

The notation: F,=d0F dp, G,=dG/dp, F,,=F/dpdq, G,,=8G/Opog; pq=—sz,y,t will 
be used. The convention that F denotes F(ryK.Ynktnax) and that G denotes G(tyKsj,YmK+s, 
tnx.) Will be adopted. The same convention will hold for derivatives of F and G. Further- 
more F),G,,F.@, will denote dF /dt, d@/dt, 


F,,F?+2F,,FG+F,@+2F,F4+2F,,4G4+ Fu, 


G,,F? t 26,14 T GG T 2G,,F T 26,4 T Ga 


respectively. 
Then it is seen that 





| ) 
F (ter. YUmkK4 stnx+-i) 
I T F Afters — Fn } T FL (Ymk4 ~“YUmkK) T F jh T SF r2(ImK4 j—Imx)* 

T 2 Fey (mk t — mK (Ymk tj YmkK ,y Fy (Ymk t — Ymk )? T Fgh )° 

+ 2F i (tme+j—Amx) (Gh) +2F y(Ymasi— Ymx) (Gh) | 4 
G(tmKksY mK tmK) = G— G(tmx4j—Fmx) — Gy (YmK+j— Ymk) — G(jgh)+ ... 
= — pin via os — , (Ah)* - ai (3.1) 
S,= A, KAF + (Kh)0, [FF + F,G(txmYnnstmx) +F lt (An) Peek 

T 2F FG (tnx Ymkstmn) T FQ? (lmx,Y mK tmK) T 2F,,F T 2F G(tmK,Ymk:tmx) 4 F | 

| (Kh °6,[F (FF st FG (tm Ynk tmx) 4 F,) } Fy (GA tmx. Ymk tmx \F 

t G (Ip KY mK tmk \G( Ink. Y mK stmK ‘7 G (img Y mks tmK ) )| T 

(n=1,4,7) 
4 
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Let y,(t) denote the solution of (1.2) which assumes the value y, at ?,. Then 
Ymk+j— Ymk = Ymk4 YmK+j(bmK) + YmK +i (tmk) — Yk: 
Since this will be a third order integration method, 


Ynksj (tmx) —Ymx=ah' + bhi 4 


Now 
er. hae — 
YmK4 Ymk+j(tn) Gjh+ “> lf (2nK+97) UmK+9) bmK-4 )+-G6,6-+-G6, 
and hence 
oe [) —— — ; 
YmK-4 Ynk=—Gjh+ - GF (tins js Umk+ js Emk+j)+G,G+G,|+ ... (3.2) 


With equation (3.2) and by repeated substitution, eqs (3.1) assume the form: 


F (amis js Ymxsjs tmxee)) =F +A KA, F+jh(F,G+ F,) 


$+ [Fe F?K?AN + Fy FGKA§ + Fie? + FPP + 2k KA G+ 2k GPL... (33) 
G(imk:Ymktmk) G— GFK Ayh ~jh(G,G@ +Gi+... (3.4) 
(Ah)° 


Y,=FKA,A+ (KA)?T, Fi + = 4.f 24 (Kh)*¢, (FF, + F,G)) 


~(Kh)?T,|GFK AA + G,Gjph+Giyh\> ... (n=1,4,7). (8.5) 


4. The First Approach 
The difference YmK+j(bmK+j41) — Ymk+541 will be expanded in a Taylor’s series and the coef- 
ficients of all terms of third order or less in A will be equated to zero. The resulting equations 
will be used to determine the integration and extrapolation parameters. Now 
YmK+j (mx j41) — Yk +541 = Ymxs (lms j41) — YmK+ A al, — ad atl, 


With eqs (8.3), (3.4), and (3.5) it is seen that: 


y ] » y , ’ y ’ 
YmK+j(bmK+jo1) = YmKk+j7t Gh SIGE (ting Ynk+jstmK+j) + G,4+ G,| 


l ‘ — — . sn 
; glelGn t 26,6 t GC? t GF? (tar. YmK- eink ) 
) 
t+ 2G6,F (ams 1s Yk +jstmn+j) + 2G ry GF (tmx Ynk+jstmK-4 
l ’ ’ ’ y y U 
t gelG (GF (rms YmK+4 stink \+-G GG, 
, > (4.1) 
- Gf F’.(2mx-4 Ymk+ itm) FP (tins jYmK+jtmn 


T GF, (tmx: YmK+ stmk-4 7| PF Alums Ymk stimk 


’ l » y l ’ 
YmK+ 571 Gho (VG, t gl Ge 





T gli(GG, T (fF) T oh’ G,(f IK A,-+ FG) +- ft fe oe 
) é 
do(}) AG. (4.2) 
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di()) AG h(G.d, t G,Guyh t G ypyh) 7 
1 ph|G,,25 T 267,,GZ gush T 267, X guyh T GG? (uh )°+ 26), ,G(sayh + Gy (uyh)*| 


hG4 h? (GFE Ag 1 G Gu, | Gr if) ) > (4.3) 
+ 3h8(G,,F?R? AR +-26,,GFK Agu, + 26,,FK Agu, +G,,Cu3 





+ 2G, Gui + G, 3) +hReCYG,F, +... J 
ds(j) =hG+h(G,2,+G,Gush+ Gunh) 
+ Bh[G,.32+2G,,GZ sph +2G, Zeus + G,,G2( ph)? +2G, (poh)? + Gei(poh)?) 
L ush?G,(G,2.+G,Guh+ Gh) 4 
hG+h1G,FK A;+G,Gu.+ Guo] 
3h ep F 2K? A2-+20,,GEK Nous +2Ge,FK Napty + Gud +20, pd + G rand] 
| uh®G.(G,FK Apt G,Gu,+Gyy)+WACTIG.F +... (4.4) 


The following system of equations results when the coefficients of the various terms are 


set equal to zero: 


» 
(4.5) AQ) + ay + a.— 1 Oty Mi + toms 

l l 

2141 7 A2he— 5 OM Mg 6 
, . tie ‘ ] ’ 
(4.6) mA K+ aAK = cpg ACA, 6 (4.9) 

2 ) 
: ite aetna ae — | 
(4.7) ok? an Az kK? 2s a, Kk? 477 ak?! 7= “a t5 (4.10) 
} > 2 
. : i - 1, KA 
(4.8) 0 AV Ag+ cop XK Ay=-~ oe KeAT+-aeKe-r, =@ — (4.11) 
} ) 2 


Equations (4.5) are the usual equations for the integration parameters of a third order 
Runge-Kutta method. From (4.9) it follows that A Ay=. and hence by (4.6) KA;=yp». 
From (4.10) and (4.11) it is seen that AA,=j. Equations (4.6) through (4.11) may then 


be replaced by 


My Mm . ) Rod le , 
A, A, K aT y+aQl; jae (I +380) (4.12) 


No attempt will be made to discuss all possible solutions of (4.5) and (4.12). However, 
a few obvious facts will be pointed out. From (4.5) it is seen that », #0, u:40, a0.  Like- 
wise y, #0. It is clear that A,, Ag, A;, 'y and T; are independent. Since the matrix 


is of rank three any set of third order Runge-Kutta integration parameters may be used. 
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One could at this point make an analysis of various types of simple solutions of (4.12), 


However, in the next section it will be seen that there are some other considerations. Two 





sets of parameters are given for comparison purposes. ) 
"pa ee - P ~\\°) 
Qi, Oy, Oe, My, Me, My (arbitrary solutions of (4.5)) 
j My Mo 
te K M K 7 kK Ne > (4.13) 
1+37 \ 
= — TS Ng =Ag—=Ap=Ag—=Ayp=0 
bay, AK } 
0 3 l 7 
fa i a | ck 4 
2 
)1=-; pw, = 0 My | 
») 
> (4.14) 
} 2 2 : 
A, ==> Xr: -—X d- ie 
2 ' 3A Vy, A 3 
Ae Ae==Ag== As Ag—A_— 0 7 





5. The Truncation Error 


In this section the truncated terms of the various expansions will be considered in some 
detail. Atermshall be said to be 7(j* A®%h’) if it is of the form 72 A°h?/ where / is a function in- 
dependent of 7, A and hk. It would be desirable for all of the fourth order terms truneated in | 
the integration of (1.2) to be T(h*). The procedure would be pointless if some of the truncated 
terms are T(A‘h*). It will be seen that there are terms which are TUvh*) and 7(7°A7 ‘A ) which 
cannot be simultaneously eliminated by any choice of extrapolation parameters. 

All of the truncated terms from yx 4) (teas jc1) —Ymxsjat Which are not 7(h*) are listed below. 

It is assumed that A;=j/K, Ay=w,/A, Ay=mw/K. The terms from Yrs j(tewsjoid, Ci(y) and 


d,(j) are grouped in that order. 


NE i es ET, .. « 
4 (G,,F + G,,G t G,,) by + : D4 G,F 
(1447) , 7 gr , (14+4j+67) , » , (1+4j—-67) 1 oa. 
+ D4 GG, 1 = 2 D4 Gf ,+ 24 Gf G43 


w ACT, (G6,,F + G,,G+G6,,) P+ K°o,G,F, Fy + bK Ay G, fb, + K?(Koy— jl) G, FG); 
be k? 3(¢,,F t G,,G T GF; t K°¢o,G,F, F 1 1K A,G, F, 1 k* (Ko; } rs) G, BG, a AC,G6,6 IK. 


The original statement was that y(f) did not depend strongly on s(t) or varied more rapidly 
than s(t). This statement may be replaced by the following explicit assumptions on / and G: 
Assumption 1: The following inequalities are valid: 
F’,| <|G,|/K, p=z,y,1,2. 
This assumption implies that y(f) does vary more rapidly than «(¢). 
Assumption 2: G,is T(K~"'). 
It is seen from (4.12) that Ty and Tz are TyA~* 
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With assumption | and the preceding remark it is seen that there are two groups of terms 
which are not 7(h*). They are 


1 +-4)+-6)*) , l prs lie 
- G.F.,, 5k ‘A,G.F>, 5k ‘AGT; 

and 
(1+-4)7—6)") 
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G, EG, Ke Ko, IVDGF YG, . k*| Ko; —jl; IGF G,. 


These terms would be eliminated if the following equations held: 


1+4j+6)? 


12h aA, T atyA; 


1--4)—6/° 


24K? a (Ad,— JT y) +e2(Ko;—T)). 


These two equations may be transformed into 


1+4)+-67 
oR Z aA, T QoAy;. (5.1) 
1+8)+6/° 59 
4K =O QD, t QoQ,. (9.2) 


Unfortunately, it is not possible to satisfy (4.12), (5.1) and (5.2) simultaneously. The 


equations 


: , 143) 
aT sy+ aT; ar 76) 
1+4)+6/? 
ay As T QA; = Te : 
1+8)+6/° 
194 AD, = WK? 


are incompatible. It is possible to satisfy (4.12) and either (5.1) or (5.2). Since the terms 
involved by (5.1) and (5.2) are T(PAT A) and T(AK'), respectively, the smallest truncation error 
results when (4.12) and (5.2) are satisfied. 

If (5.2) or (5.1) are to be satisfied, then some of the simplicity of the extrapolation coeffi- 
cients is lost. Two sets of parameters are given below; the first satisfies (5.1) and the 


second satisfies (5.2). 





Mi, Me, Ys Yo. Ys (arbitrary Runge-Kutta integration parameters) ) 
) My 1+ 3) AsY2 
A — >—A-—Ag \;=- > — 
hk a 6ayy,K? 1 . (5.3) 
1+4)+6/7?—2y7,A(1+-3)) 
ha = = Ng === 0 
12Kayy2(¥2—71) K J 
i, Moy Yi. Yo. Ys (arbitrary Runge-Kutta integration parameters) ) 
] ; 1+3) = eyo 
d, J r, ~ a e , \;=-= = Y2 
K K sayy? ¥, (5.4) 
1+8)+6/7 ' 
¥ . “J | he at he A, Ae a 0) | 
24k ‘aan yy kK 











These equations are somewhat more complicated than those given in section 4. 
With assumption 2 it is seen that there are three groups of terms which are not 7(h*), 
They are 
(1+8)+-12K"9T;) , nr iets ater eanal ee 
-—" YGFF,, K°OGFF,,  K°¢G.FF, 
and the two groups found with assumption 1. These terms would be eliminated if the following 
equations held. 
14+4j+6;? 
4K : a As T OtA7 } 
1+8)+67? 14-8)4+12AK°P, 
<= <—_s “=ay4+ ayy. 
24K" 24K" on 
Hence if T, is set equal j?/2A® the same equations are found as with assumption 1. 
6. The Second Approach 
In this section the extrapolation parameters will be determined so as to make the trunca- 
tion error of the extrapolation as small as possible. This procedure is similar to the analysis 
for the Runge-Kutta integration of one equation. 
The extrapolation of x(t) to any value ¢,,< to any value f,,.<+-7 Is given by | 
mK tT No (7) ko T i, (7) hk, T ho (rt) he 
where ky, k, and ky are from (2.1). Ap(r), Ay(7) and Ay(r) are determined by equating the co- i 
efficients of A in the expansion of the error equal to zero. ~The resulting equations are 
Ay (7) +A, (7) +AQ (7) =7/Kh (6.1) ; 
D  giegedeneee : ae 
M(t) yi t+ Ae (7) Y2=5 77/( KA)? (6.2) 
‘ ow 5 pcan . 
di (7) v7 +-Ao (7) ¥3== 7°/ Kh) (6.3) ) 
”) 
l swe , 
ho (tT) V3 g 7 /(Kh)’. (6.4) 
) 
Equation (6.1) results from equating the coefficients of A to zero, (6.2) results from equating 


the coefficients of h? to zero, and (6.3) and (6.4) result from equating the coefficients of h* to 
zero. In general only three of these equations can be satisfied One would naturally choose a 
solution which satisfies both (6.1) and (6.2). 

These equations may be written for the extrapolation parameters. The resulting equations 


are: 
/ My Me as 
A = A > A; > (6.5) 
A , K K 
~ ( os \- ( ‘| )? 
I, Te r, mc —Tr;, lr; “aR? —T, (6.6) 
y° (7+ 4,)° (j}+- po) _ 
A — A : =——A A; : =, A, (6.7) 
' 3K . 3K ; 3K? a 
lh (j+m)) (j+- m2)" . 
=—5 : = — 7 : = —?,. (6.S) 
= 6K? me bk °1 %7 6k % 
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Equations (6.6) and 6.5) imply that (4.12) is satisfied. The converse is not true. It is easily 
shown that (6.7) implies 


6)? +4) + 4 (ami + apy) 


a Ay T QA; 12h3 


and that (6.8) implies 


6)? J 4) } 4(a,p} | Qt) 


Py Aa; 4K 


Therefore (6.7) implies that the G,F, term is 7(/?A~'h*) and (6.8) implies that the G,FG, term is 
TUKh*). Hence (6.7) and (6.8) have the same effect on the truncation error as (5.1) and 
(5.2) although (6.7) and (6.8) do not actually imply (5.1) and (5.2), respectively. 

It is rather surprising that [';—j*?/2A* does not appear among the equations derived in sec- 
tions 4 and 5 with assumption 1. It is certainly plausible to include this equation in any set 
of equations taken to determine the extrapolation parameters. 

Two sets of parameters are given, the first satisfies (6.5), (6.6), and (6.7) and the second 
satisfies (6.5), (6.6), and (6.8). 








My, Mo, Yi, Yo, ¥3 (arbitrary Runge-Kutta integration parameters) . 
27)—37.K}? 2p—3KP 
A, z—Ng —X3, A» , . 1d IF - mt 
rN 6K x (11-72) 6K y(¥2—%N) 
ree eee oot | aoliiie dla 6.9) 
 .. nee r. 2(j+ m)°'—3y2K (7-4 my) ‘ 2(j+m)°—3K (j+m)% > ( 
—_ 6K y,(11—Y2) : 6K *y2(¥2—71) 
»— M2» ‘ 2()+ m2)*— 32K (J+ me)" _ . 2(j+m:)°>—3K(j+m)?1 
 — . 644 y,(¥1—2) = 6A y2(%2—71) J 
My, Mo, Vi, Yo. Ys (arbitrary Runge-Kutta integration parameters) 
} ? YJ" 7 
r ~>—As—Az, Ay heey oe he ® Za 
= oh *y, 6K 27;  6IONYs 
(3 ae ees 3 (ij 3 7 (6.10) 
My p> Bi)” I+ Bi J+ Bi} 
A > —As—Ag, A; “de, 62 — . 
at 2h*y, HA yy; . 6K ny; 
Me (j+-m2)? — Y¥2(J+M2)* (J+ H2)* 
Ay =F — Ns — Ag, = ete 7g a 
AK 2h ae 6A 'YIY3 P 6A “Vis e 


If only (6.5) and (6.6) are to be satisfied there is considerable simplification in the extrapola- 


tion coefficients. The following set satisfies (6.5) and (6.6): 





r J —dX nN J nN 0 | 
, eo * 2K*y, , 
My (j+m,)° 
nN -—A; A; —__— —X> A,—0 (6.11) 
a ° 2k*y, : ; 
Me (J+ M2)? 
ry ~—Ag Ag 5 -\. A,— 0 
‘ Kk 2K *y, . , 
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7. Fourth Order Formulas 


An analysis similar to that of sections 4, 5, and 6 may be made for a fourth order integration 
method. Only the results of such an analysis are given here. The equations for a fourth 


order integration method are: 


ko= KhF (tmx, Ymk, tmx) 

ky = KhF (tan +rniko, Yax t+ yiho, tnx t+ Wh) 

ko = KRF (tinge t+¥3hi + (Yo—-Y3) Ko. Ym + Yahi t+ (Y¥2—-Ys)ho, tmx + Y2kth) 

k3= KAF (tne t+ok2+ski + (¥s—¥s— 8) ho Yk + Yoh2t+shi t+ (va ¥s— Vahey tme +skth) 
hiyp= KG (funk, Unks tne) 

hi, =KhG (Qn +rniko, Unk +Miho, tax + yh) 

ho=KhG (tan t+ryakit (v¥2—Y3) ko, Yak + shi + (V2—-Ys)ho, tmx + ¥2kth) 

hs3= KhG (ame +yokot+yskit (¥s—¥s— Yo) ho, Ymx t+ Yolt2t+ shit (v¥a— 5 — Yoho, tra velth) 


DL ¢m4+1) K — Lmxk Tt Bok'o T Bk, T Bok, T Bk, 


4 
mK + —>> ACD; 177 mK 


J — 
i=] 


dy(j) WG (tiers Ymk+jtmK+;) 


~ 
di(j) AG (Imesjt> i,k, 594 mK+j71 py ,t, K+ij71 uh) 
i=5 


dy(j)=hG (tmsjA DS Aki—9, Ym s+ Mali + (Me Ms) astm Meh) 
i=9 
; eo 
d3(j)=hG@(amnsj+ > Ak, 133YmK+j 7 polly T asl T (Ma Ms Me Lo stmns T gh 
i=13 


Yk +541 = Ymw+j) + coll + a4, + aly + ads 


The integration parameters must satisfy: 





Ay Oty + Ay Az l Cy My t Ot» iby t Qn Ms , 7 

1 1 
Ot) My TT O2M2 TT A3h4— 2 Of) My Mo M3 -t- Og Mo M4 Met Og My Magis 5 

> 7.1 
2 2 2 2 2 2 
Ot) MT ln Me T Og My 3 O2Mi Ms TT O3MoMa TT O3éM1 Ms . 
1 

Oly hy Mgt Oly Mo Me + Oly My Ms 6 C63 My MyM = 94 J 
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The extrapolation parameters must satisfy the following equations in order for the 


integration method to be fourth order: 





j : , 143) 7 
A, "Kk a;4 al +- a3 1" 3== 6K? 
im 1+ 4) { 67° 
A. Ke a, A; t Oty Ay t OyA)3 lok 
Me : . r 3+87 — 
Mn= Oy My P54 opto Py 4- guy The aa (7.2) 
M2 1+-8)+-6/? 
Ais “Kk O)Ps5 + AsDy + AgH 3 4K 
es ; ’ , 144) 
] l “9K? Qops | a7 Ops | s—T Aap | ) ~O4hK? 7 
where 
A ” A, l Nn yi T A, +271 Ans 3 A, viv, +1 T V5An+2 T Virn+3 
r’, ViAns T V2An+2 T V4Ans 3 Py ViVsA,. 42 T (NYS 275) Ans 3 


The truncated terms have not been examined as in section 5 due to the extremely tedious nature 


of such an examination. 
The equations analogous to (6.5), (6.6), (6.7), and (6.8) are: 


A, -Z he x Ay ‘R A= (7.3) 
— . (j+m)? (j+m)* ,, (j+us)? 

"3K no eae 7 

(74 ) (7+ yn.) (7-4 )3 

aie) i) ce) cr) ire va 

] { ) + My) (j + My) (74 by)° — 

ws ta) +e)? +H)” _ 6) 
*16K3 Ps bk ° ? 6A ?1 P13 6k? % G 


The relation between the A, and the A,, T,, A, and @, is 


\ 1 | | | r, 

I’, oO ¥ Y2 Ys Anat 
A, Ov vi An+2 
>, O O ys (Ns V2¥68)) | Ants 


This relationship can be inverted and 


A, \ 

An+1 | 
A 

An+2 A, 

= ra) 
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where A is the matrix with elements a,,/p, 
p= Vil (¥2— M1) ¥2(¥1 ¥5 + V2 Yo) +11 ¥3 Va(M1 vs)] 


aj,=p yp (V1 = Yo) (11 V5 V2 Vo) — V1 V3 (¥1— Vs) = 3 = (V2 — V1) (M1. V5 V2 Vo) V1 Ya (V1 — ¥4) 


14= (V1 — Y2) (V2 Ya (M1 Ys) 
dy, =0 22 = YolVi Vs + Y2 Ve) — V1 3 ¥43 23 = ¥1 V3 Vs — Y2AN1 Ys + V2 Yo) 23> V2V4 Ve's 
3, =0 (132 Vi(¥1 ¥s +7270) 33 = V1\.'71 Ys 1 Y2 Ye ) 33= 71 V4(¥1 — 94) 
3 2 ~y. «4 
ay =0 (42 Y1¥3 43 71¥3 (144 ¥1 ¥2(%2 Yi) 
If Runge-Kutta-Gill integration parameters are used, then 
(| 3 2 0 7 
0) 0.58578646 2.8284272 6.8284272 
A 
0 3.414214 6.8284272 6.8284272 








0) —1 2 0 ; 


8. Experiments 


Experiments have been made with these formulas on three pairs of differential equations. 


They are 


4 


dr 
dt 


x/2 r(Q0) | 


(S.1) 
dy ~~ 
J — x eos (251) y¥(0)= 1/1250. 
di rco ) Y Wd 
dr 
dt 


=y° r(O0)=a®/13 
(S.2) 
dy lou (0) =a 
: } ( 
dt a(t+1)” 4 
dr 


.1+ (1.5 107~*) cos (1507) r(0)—0 
dt 
dy 


rcos (257) y(0)— .OOO16. 
dt 4 





4 


Equation (8.1) is a typical example of a system suited for split Runge-Kutta. Equation 
(8.2) is an example of a system not suited for split Runge-Kutta. Yet even in this situation 
there may be some advantage to using it if the solution for y(t) is needed with much more 
accuracy than that for s(t). In (8.3), y(t) is much more rapidly varying than (ft), vet the 
small oscillation in x(t) gives it a comparatively large second derivative. The experiments for 
each of these equations are discussed in considerably more detail in sections 9, 10, and 11, 

No experiments were made with parameters derived from assumption 2 in section 5. 
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9. First Example—The Typical Case 


The solution of (8.1) is 
r(t)=e'? (9.1) 
y(t) = (4 cos 25t+-25 sin 25t)e"?/625.25. (9.2) 
There are three basic variables in the solution of (8.1). They are (i) the method of solution, 
(ii) A, (ii) A. The methods of solution are numbered 1 through 6 corresponding to equations 
(4.13), (5.3), (5.4), (6.11), (6.9), and (6.10), respectively. The same integration parameters 


were used throughout. They were: 


a By 4 q. 7) My ] z. V2 Me i | 4, 3 Ms 3 4. 


~~ 


ay Bo 2 9, a, By ] 


In figure 1 the errors in the y integration are studied as a function of K for various methods 
of solution. The integration interval / is fixed at .01 and the error is the maximum error for 
O0<t<1. In general this occurs close to t=1. The value of the error of the z integration 
at f=1 is also given. Here the variation of A is equivalent to varying the integration step. 
This error behaves like A®? in the range of values considered. 
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There is very little difference m the results for the first three methods and only one curve 
is given. These methods proved to be inferior by a large factor compared with the second 
three methods. 

There is little difference in the results for the fourth and fifth methods. Only the results 
for the fourth method are plotted. This lack of difference can be explamed. Method 4 
6.2), and method 5 corresponds to satisfving (6.1), 


corresponds to satisfying eqs (6.1) and 
(6.2), and (6.3). Now (6.3) arises from setting the coefficient of Fy equal to zero and for this 
system F=ar and F,=0. Hence there is no improvement when method 5 is used rather 
than method 4. 

On the other hand method 6, which satisfies (6.1), (6.2), and (6.4) is a definite improvement 
over method 4. Equation (6.4) arises from setting the coefficient of FF, + F,G,=e'?)8 equal 
to zero and hence an improvement would be expected. 

There is another phenomenon apparent in figure 1 and a very significant one. That is the 
fact that the error does not decrease with A after a certain point. Indeed, the same accuracy 
for y(t) may be obtained with Ah=0.35 as with Ah=0.01 if method 6 is used. This phenom- 
enon Is explained by the fact that there are two factors that contribute to the error in the 
y(t) integration. One of them is the error in extrapolation of (ft) and the other is the error 
inherent in the actual integration of y(t). For some values of A the extrapolation errors are 
dominant and for other values the integration errors are dominant. When the extrapolation 
errors are no longer significant then improvement can be made only by decreasing h. 

In figure 2 the errors are studied as a function of A with A fixed. This was done only for 


FiGgurReE 2. Marimum integration errors as 


a function of h and the extrapolation 
method, 


Wy is fined at 25 
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the first and fourth methods. 


Figure 


») 


is self-explanatory except to point out that things are 
not always as one thinks they should be in numerical analysis. 


There is no particular explana- 


tion for the fact that the error does not decrease monotonically with h, it just happens that way. 


Figures 3 through 7 are plots of error curves for various situations. 
Figure 4 is for method 4, A 
It is identical with the curves with A=5, 1. 


Kk 
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Ol and it is identical with the curves with A 
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K=35. Figure 7 is for method 1, A= 25, h 
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Figure 3 is for method 1, 
Figure 5 is for method 4, 
O- 


a, 


OL. 
Figure 6 is for method 6, K 
10, 5, 1 and very similar to the curve with 
1/W2 = .0084. 


25 and h 
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Figure 6. Error curve for method 6 with 
K= 25 and h=0.01, 


Identical curves are obtaired for smaller values of K 
ind for A’=35 the error curve is very similar 
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Figure 7. Error for method 1 with 
K= 25 and h =0.0084. 
It is seen that the relatively larger error here is due to 
the chance shape of the curve and not to any basic lower 
ecuracy. 
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10. Second Example-—-An Unprofitable Case 


The solution of (8.2) is 
r(t)=a*(t+ 1/13 (10.1) 


y(t) =a(t+1)?. (10.2) 


All of the results discussed below are for a=.1. 
In figure 8 the errors are studied as a function of A’ for various methods of solution, 
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INCORRECT, ; oa 
METHOD6 y / Pl / 
/ \ Y a 4 / 
f ie. P 
i he or 
Pd 
a J ‘ f= METHOD 4 +5 
/ ~ 
/ 4 ~ : 
a rf / “~ METHOD 6 o 
r~ ~k / ¥ / \NCORRECT 
/ ae METHOD 6 
In XL / 
/ i / Figure 8. Integration error at 1 0.2 asa 
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Methods 4 and 5 give very similar results while method 6 
Also included is an “incorrect method 6” which 


In this method Ag was computed from 


Methods 1, 2, and 3 are not used. 
does somethat better than either of them. 
gives the best results of all for larger A’s. 


() T M2)” y2() + My)” 
s 2h*y, OA viv: 


instead of its true value given in (6.10). 

The fact that method 6 is better than methods 4 and 5 cannot be explained as previously. 
Indeed, by that argument method 5 should be the best and method 4 the worst with method 6 
in between. Reeall that method 5 satisfies (6.3) which is obtained from setting the coefficient 
of F, equal to zero, likewise method 6 is obtained from setting the coefficient of FF, t FG, 
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ERROR 


ERROR 


Now in this case F,=120a* (t+-1)" and F,F,+ F,G,=6a® (t+1)" and hence 


equal to zero, 
This points out that the truncation error is not always a good 


method 5 should be the best. 
measure of the integration error. 

It is seen that the integration error does not level off as AK decreases. 
fact that the major portion of the error in the y(t) integration is due to the errors in the inte- 
gration of x(t), as differentiated from the errors in the extrapolation of x(t). This is borne out 
In figure 13 the y(t) integration error is plotted for the system of equations 


This is due to the 


by figures 9 and 13. 


dz =@°(¢-+-1)!? 
a7. 

dy -12~ 5 12 
dit -l3r/a°(t+-1) 


which have the same solution as (8.2). Method 6, K=20 and h—.002 is used. This leads to 


much more accurate values of 2(t) which in turn lead to a much more accurate y(t) integration 





ERROR 


ERROR 


even though the extrapolation accuracy has not improved. 


In figure 9 the error curves for methods 4, 5, and 6 


and for 0<t<.2. 


The error curve for method 4 with AK 


Lh 


are compared with K=20, h=0.002 


0.002 is given in figure 10. 


In figure 11 two error curves are given which have the same interval of integration for 





















































the «(t) equation. Both use method 6, one with A=1 and h=.02 and one with A=10 and 
| ore: | 
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Figure 9. Error curves for methods i © ee 
and 6 with K= 20 and h=0.002. 2X10 > . 
The improvement of method 6 seems to lie in the fact , 7 
that the error remained positive much longer than for Figure 11. Error curves for method 6: with 
the other t ’ ethods - : ; ’ 
le other two meth K and h taken to be (10, 0.002) and 
(1, O.O2). 
These curves emphasize the fact that the accuracy is 
determined by Ah and not h alone 
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Figure 12. The unusual error curve jor the 
Figure 10. Error curve for method 4 with incorrect method 6 with K=10 and h 
K=1 and h=6.002. 0.002 


rhis corresponds to the norma! Runge-Kutta integra 
tion procedure 


The improved answers are apparently due to the fact 
that the error remains positive for a long period 
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Figure 13. Error curve of the modified 


> equation for method 6 with Kk 20) and 
h—0.002 
- : 
/ \ Che modification removed the dependence of the 
L st 4 4 / ‘ equation on yi Ther 
if i) 


factor ¢ 








h=.002. The close similarity of the two curves again emphasize that the accuracy of the 
r(t) integration determines the accuracy of the y(f) integration. 

The unusual error curve for the “incorrect method 6” is given in figure 12 with A=10 
and h=0.002. 

Some studies were also made with a=2 but the overall behavior was the same. 


11. Third Example Noise 


The solution to (8.3) 1s 


r(t)=t/10—+ 107° sin 150 


ty(t) =4 107 [sin 25t—6 sin’ 25t— (48 5) sin’ 25t— (32.7) sin? 25t! + 0.004 |f sin 25t+ .04 cos 25¢]), 
This system Is vers similar to (S.1) with the addition of a rapid oscillation on «lt However 
the experimental results bear no resemblance to those for (8.1). In facet they seem to be com- 


pletely chaotic. Although the errors in the y(f) integration are all relatively small and fairly 
uniform in size, there appears to be no correlation between the maximum error and A’ or method. 
An examination of the error curves reveals the same lack of trend. Even the errors in the 
r(t) integration do not behave rationally as a function of AY which is equivalent to the inte- 
gration interval. Indeed the largest error occurs for A= 4. 

In figure 14 the y(f) integration errors are plotted as functions of AU and method. The 
v(t) integration errors, which are independent of the extrapolation method, are also plotted 
as a function of A. 

Figures 15 to 18 are error curves for various situations. 


12. Numerical and Efficiency Considerations 


The ingredients of an ideal situation for these formulas are as follows: (1 Fir.yt) is com- 
plicated and difficult to evaluate, (2) the solution for y(f) is relatively insensitive to errors in 
r(t); ie., the main error source is the inherent inaccuracy of the y(t) integration. The first 
ingredient magnifies the efficieney and the second is necessary to attain high accuracy. The 
second ingredient may be replaced by “the accuracy desired for y(t) is lower than that required 
for r(t).”’ It should be remembered that numerical integration is a delicate business in general 
and these formulas are particularly so and hence each case should be examined on its own 
merits. 

A comparison will now be made between the work required for split Runge-Kutta and 
the normal Runge-Kutta methods. 

Unless the parameter AK is changed very often the extrapolation parameters should be 
computed only once. For the integration of (1.1) and (1.2) from ty to ty + Ah the following 
tables give the required number of operations for normal Runge-Kutta and split Runge-Kutta 
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Figure 17. Error curve for method , with 


K=5 and h=0.01. 
The error oscillates as in figure 16 except it is now 


scillating about the line; error 1.2/+-0.2)10-°, instead 
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methods. Table 1 is for third order integration and table 2 is for fourth order integration, 
From these tables it is seen that if a large A may be used then there is a large saving in com- 
putation even if Fis extremely simple. If Fis complicated then the efficiency is even greater, 


PABLE | TABLE 2 
Add Mult F eval (i eval- Add Mult F eval G eval 
uation uation uation uation | 
Normal Runge- 
Kutta ; 2h 2A 3h 3h Normal Runge- } 
Kutta SA WA ih ik 
Split Runge-Kutta 
method 4 14(A+1 1TA+1% 3 3h+2 Split Runge 
Kutta 2A +27 20h +34 1 1Ah+3 
Split Runge-Kutta 
methods 5, 6 16A+14 In A+1 3 3h+2 
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A Reduction Formula for Partitioned Matrices 


; (April 


A theorem of L. Goddard and H. Schneider, concerning square matrices 


orders n and m, respectively, which satisfy ar 
is generalized here 
which satisfy AN,—N,B, where X; has dime 
find reduction formulas for partitioned matri 
n,n;, and satisfying equations A,;N;— NX, 
also generalizations of a theorem by J. Will 
submatrices are all square and satisfy AN 


1. Introduction 


| Given a matrix A=(a,,) of order n, 
,a nonsingular matrix ? such that 
; 


BOO 


A=P-'AP (1) 


'where Bis a square matrix of order r, and @ is an 
rX(n—r) matrix of zeros, A is called a reducible 
matrix and the formula (1) is a reduction formula. 
The characteristic equation of A is then factorable: 


-- 


A—Nl, B—xJ,||\C—Xx, 

where l, represents the identity matrix of order n. 
The results in this paper concern a reduction 
formula of the type (1) for partitioned matrices, and 
connection between, and an extension of, 
Goddard and H. Schneider,’ and J. 


vive a 
results by L. 
Williamson.? 


is generalized in theorem 1, and this is applied in 


partitioned matrices, 


A, A, 1, 
] bes Ass L., (9) 
Ayn Ag ... Ag 


a 


+ A! ee 2 
Theorem 2 contains the results of Goddard and 


-__—— 


| L.S. Goddard and H. Schneider, Matrices with a nonzero commutator, Proc 
"amb. Phil. Soe. 51, 551 (1955 


J. Williamson, The latent roots of a matrix of special type, Bul. Am. Math. 


Soe. 33, 585 (193) 





for rectangular matrices A 


if there exists | 


The theorem of Goddard and Sehneider 


theorem 2 to give a general reduction formula for 


where the submatrices .4,; have dimensions n,X7,, 


Schneider and those of Williamson as special cases. 


Emilie V. Haynsworth 


13, 1960) 


A and B, of 
1 equation AX = NB for some nm matrix X, 
and B, with dimensions n,X no, m,X mo, 
isions n,m, fori=1,2. This result is used to 
ces With submatrices, A,;, having dimensions 
B,,. The reduction formulas given here are 
iamson concerning partitioned matrices whose 
XB, where B is triangular and X is square. 


2. Results of Goddard and Schneider 


Goddard and Schneider showed that if two square 
matrices A and B of orders n and m, respectively, 
are related by the equation, 


AX=XB, (3) 


where Y is an 7 m matrix of rank r, then there 
exist matrices P and Q, depending on .X, such that 


E+ k O 
P-*AP , QV ' BY , (4) 
oe ¢ * DPD, 


where / is square, of order 7, so that A and B have 
r roots in common. 

More generally, Goddard and Schneider proved, 
for any m*n matrix A and any polynomial f (7, y), 


S( EG) * 
P-'f(A,XK)P 
0) FCO) 
(5) 
f( kG) O 
OU BANG 
* f(D,0) 


where G depends on K. So all such functions of A 
and B are reducible if r<min (m,n), and each cor- 
responding pair (for which f and A are the same) 
has r roots in common. 

If r=m<n, A and f (ANAK) are reducible and 
(1) holds for A and B as defined in (3) and C as 
defined in (4). 

If r=m=n, for the matrices in (3), we have 
N'AN=B, so if B is reduced, A and f (A, YAK) are 
reducible. Naturally, similar results would hold 
for Bif r=nSm. 

Thus, although the results of Goddard and 
Schneider concern a pair of matrices, the cases given 
above can be considered as leading to a reduction 
formula for one or both of the matrices. 


171 








3. Results of Williamson 


J. Williamson? deals with partitioned matrices 
in (2) in which all submatrices, or blocks, are 
square. He showed that if a partitioned matrix A, 
of order nt. has blocks of order n which can be 
simultaneously reduced to triangular form, then A 
is reducible to a block -triangular me LtriX, Le. a par- 
titioned matrix (2) in which A,,—0, 7 >i. 

Williamson’s results could be expresse “i as follows: 

Given a partitioned matrix A=(A,,) of order nf, 
with blocks of order n, if there exists a nonsingular, 
nn matrix VY such that 


as 


AA, ,A=B 
or 
AyA=XB Gj soe oe t) (6) 
where B,, is triangular with elements A, . 2 2, 4 
on the diagonal, then A is similar to a bloek-trian- 


gular matrix with blocks 


on the diagonal. 

This theorem is then generalized by Williamson to 
show that, given any partitioned matrix A= (A 
satisfving (6) with B,, triangular, the partitioned 
matrix which has blocks 


G j f (A, ,;) 
where /,,(A) is a rational function of A with non- 
singular denominator, has as roots the roots of the 
n matrices of order ¢, 
Gy =f (0s (k=1, ,n; t,j=1,.. .,f). 


4. Connection Between Theorems of Goddard 
and Schneider and of Williamson 


It can be seen, from eqs (3) and (6), that although 
the theorem by Williamson and that of Goddard and 


Schneider are quite different, there is a connection 
between them. 
In theorem 1 we generalize the results of Goddard 


and Schneider to rectangular matrices ol and B of 
dimensions n,n. and m,<m,, respectively, satis- 
fying 

AX,=X,B, (7) 


where YY, is n,m, of rank r,, ¢—1 or 2. 
We then apply this result in theorem 2 to square 
partitioned matrices A and B, with rectangular blocks 


ny Xn; and m,*m,, respectively, in which the fol- 
lowing relations hold between the blocks, 
AulAs=AwBu, 3=1,...,! (S 
where YY, is n,m, of rank r,. 
?J .Williamson The latent roots of a matrix of special type, Bul, A Mat 
Soc, 37, 585 (1931 
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Then the theorem of Goddard and Schneider can| * 
be obtained as a special case of theorem 2 when w, | 
set 7 : in (8S). If we have 7 

_and B,, triangular, 7, j=1, . 


the Fhe bes “mm of Williamson 


Ut Vi Ni, = | 
f, we hay 


5. Theorem of Goddard and Schneider fo! 


Rectangular Matrices 


The proof of theorem 1 is essentially the same ais} 
that of Goddard and Schneider, the only differency 
being the addition of subscripts to correspond witl} 
the two different matrices, .Y, and LY. 

THeorem 1. Jf A and B are rectangular matrice: 
satisfying (7), there he matrices P,, Ps, ©, and ¢ 


such that 
kK * | 
Ge €) 


MI! 


(4 
‘ko OO 
W ‘BQ, ( , 
* DD 
where Io is an ry ry matris, and Cand D are rectan: 
gular matrices having dimensions (iy ry) * (e—r) 
and (mm, —r1) (m:—"s), respectively. Moreover. 
Mr, yy) wa polynomial in two noncommutative indete } 
minates (which is linear and homoge neous if A and i 
are not Square, ol if l= I> and K is ah arbitrary My } 
ny matrix, \ so 
fi iG) . 
PO ALX, ADP: ( , 
0) fC 0) | 
(10 
/ iG i) ; in 
, if BKN DQ. ( ) | WI 
ADO), | Ca 
or 
where G depends upon K 
Proor. For j—1 or 2. we have Y, of rank r;, so) ™ 
there exist corresponding nonsingular matrices, P| gl 
and Q,, of order n, and m,, respectively, such that! 
Ne 
LO 
\" P iy Q) ( ) (11)) W 
0 ) | 
isan m, matrix with an identity matrix of orde! 
r, in the upper left corner. 
Let W 
mM 
~ » he iL’ o\. & 
A=P7;'AP-sz, B=0 BO,, KkK=Q;'KP, (12 
and partition each matrix so that there is an ry) > ' D 
matrix in the upper left corner, Le 1 
- « lb 
aA aA ie 
al , 
A; Ass - 








der can| 

vhen yy, 
i=] 

Ve hay, 


er for! 


sane ass 
Terene| 
idl wit! 


natrice, 
and f 


; 
rectay-| 
Na— hy 
OVE, | 
nde ter) 


and I 
{ ny 


ry, & 


°s, P| 
1 that! 
(11)) 


| 


orde! 


(12 





and Band & are partitioned similarly. 
Then, from (1] ) and (12). 


- Aj v | 
PP; 'AAY,=AYl, ’ 
A, O | 

(13) 


we 
~~ 
“ 


P'X,BQ.=Y,B 
) i) 
Since, by (7), AY:=Y,B, (13) implies 
Au=Bu, An=B»=0 (14) | 
Thus, if we let Ain E, Ae. os B, D), from 
12) and (14) we have (9). Also 
me K,, Z.. 
Pr'X, KP.=Y,K 
— 0 
(15) 
. K,, O 
QW, '(AN)QW=AY, 
ky, O 


so if we let Ay, = G, we obtain (10). 


6. A General Reduction Formula 


In this section we apply theorem 1 to partitioned 
matrices with rectangular blocks satisfying (8). We 
will in general obtain formulas of the tvpe (4) which 
ean be considered as reduction formulas for A or B 
or both. 

Specific applications of this theorem in the redue- 
tion of certain special partitioned matrices have been 
given previously by the author.’ 

Before applving theorem 1 to theorem 2 it will be 
necessary to prove the following. 

Lemma. Given a partitioned matrix of order NV, 
with n,n, blocks A,,, if 


B 0 
(02) 

* (’ 
where all matrices B,, are square, of order r, and the 
Matrices ( are (i rx (hn a> then u1 Is reducible 
as in (1) where P? is a permutation matrix, and 

B (B . €? ((' 

PROOF. The proof consists merely in defining the 
permutation matrix 77, which is equivalent to writing 
the order in which the rows ana columns of aA should 


(16) 


be arranged. 


E. Haynsworth, Applicatior 
search \ BS 63B, 73 (1959 


Let >05.4n;=N,. Then Ni=n, N,=N. If we 


arrange the rows and columns of A in the following 


| order: 
Rs 2 PF 7 
Ni, +1 N,+2 N,+7 
N,-s1+1 Ni-1+2 ~Neaitr: (17) 
/ 8 / 2 N,; 
N, / l Vv; r- 2, . \ 
N, “Tv * l, N, “se 2, ‘ N,; 


we have a new matrix A in which the matrices B, 
are Logether in the upper left corner, and the matrices 
(',, are together in the lower right corner, so A will 
have the form (1). 

THEOREM 2. Suppose we have partitioned matrices, 
A and B, with rectangular blocks, A,, and B,,, satisfying 
(8), wherer,;=r, i=1,...,t. Then, Aand B have 
tr roots in common. 

Moreore af if G 
blocks, 


(G,,, and H 


(7, ) have re ctangular 


(18) 


Giy=fhislA =f (By, Ku X;) 


where the matrices K,, are arbitrary m,n; matrices 
and the polynomials f(y) are as defined in theorem 
1, then all pairs of matrices defined in (18) have tr 
roots in common. 

Proor. Let P?, and Q, be matrices satisfying (11) 
lor tg=1, . ,tand let P? and @ be the direct sums 
of the matrices ?, and Q, respectively, Le., 


, ‘ , 
P=P,+P.4 1}. DPQ 2a°Wi. 


Then we have, using block multiplication of matrices, 


~ ~ 


P“AP=A=(A,), Q°BQ=B=(B,,), 


where 


~ 
> 
>) 


‘Ay P;, 


—_/ 
~~ 
~ 


Qi BL). 


Thus, from (9), 


—_/ 
in, 
~ _ 
~ ‘> 
~ ¥ 
eee” 

= 
* ~~ 

~ 
~~ 


DD 
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So, by the lemma, A and B are similar to the 
A _ 
matrices A and B, 
x Es ‘ k O 
A ; B 
oO € . pp 
where H=(F;,), C=(C;,), D=(D,,). Thus the 


roots of are roots of both matrices. 

Also, using (15) for the matrices A,;, P;, 2;, Qi, Q). 
we can in the same way obtain the more general 
result concerning the matrices G and /7. 

Corotiary 1. Jf a matrix A of order nt can be 
partitioned into square blocks A,,, of order n, having r 
linearly independent characteristic vectors xs, in com- 
mon, corresponding to the roots, '", h=1l,. . . . r: 
s, = 1, , t. then tr of the roots of A are roots of 
the matrices (X,""). 

Proor: Let 

Xi 


(2),J2, I,) 


be the n&r matrix whose columns are the 


vectors and 


given 


B,,=diag (Ay"”, As", +" de 


Then (8) holds for the blocks A,,, and tr of the 


174 


roots of A are roots of B=(B,,). 
arrange the rows and columns of B in the order, 


l r+] 2r+1, (t{—1)r+1 
2, rt+2, 2r+2 (f{—1)r+2 
Pie es eg acc g OR 


B will be in block-diagonal form with the matrices 
(AA) on the diagonal. 


If we now re! 


This corollary is applied by the author to certain | 


special partitioned matrices (see footnote 3). 
COROLLARY 2. /} a matrix A, of order nt. can by 


partitioned into square submatrices of order n whieh’ 


are mutually commutative and have roots, i, (h=1. 
: , n), the roots of A are the roots of the n matrices. 
Ar >, Cie he » ana o 

This corollary would follow also from Williamson's 
theorem since any set of mutually commutative 
matrices can be simultaneously triangularized. 


The author is grateful for the helpful suggestions 
of Dr. A. Ostrowski and Prof. H. Schneider. 


(Paper 64B3-33) 
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Selected Bibliography of Statistical Literature, 1930 to 
| 1957: III. Limit Theorems 


Lola S. Deming 


trices 


ertain } 





an by | 12. 1960 
; j (JanNUArV Le, FO) 
which ’ 
h- & 
trices This is the third in a series of bibliographies dealing with various specifie subjects in 
, the field of statisties. References and titles of important contributions concerning limiting 
: distributions have been taken from technical journals published throughout the world since 
ISON s 1930. 
ative 
Complete coverage is not claimed in this series of | initial page number, and the date of publication in 
bibliographies. It is believed, however, that the | parentheses comprise the reference to the original 
two prominent reviewing journals whose abstracts | article. 
. 4 serve as our source material, have selected for review Reference to the abstract: The final symbols M (for 
~ 2 +, . . . . . » ; . . ’ . > os 
Hons (the writings of major statistical importance from | Mathematical Reviews) and Z (for Zentralblatt fiir 
technical journals and publishing houses throughout | \Jathematik) followed by a volume number and a 
the world. page number refer to the abstract of the article or 
33 | This particular subject’ classification on Limit | book appearing in the reviewing journal. 


Theore Wis follows two earlier ones on Corre lation and 
Regression Theory and Time Series.' The 675 refer- 
ences have been extracted from a card file of abstracts 
taken from Zentralblatt fir Mathematik for the years 
1930 to 1939, and from Mathematical Reviews for 1940 
through 1957. This file of abstracts is maintained 
ona current basis in the NBS Statistical Engineering 
Laboratory. The abstracts are coded into categories 
of subject matter by the subject classification used 
in Mathematical Re re ws, One abstract nay be 
classified under several subjects, hence may appear 
in more than one place in this series of bibliographies. 

To transeribe the material given here from the 
abstracts, the references were punched onto SO- 


column cards thereby necessitating severe and un- 


conventional abbreviations in many cases. 
The following information is extracted directly 
from the abstracts: 


Ackermann, W. G., Eine Erweiterung des Poisson- 
schen Grenzwertsatzes und ihre Anwendung auf 
die Risikoprobleme in der Sachversicherung, Schr. 
Math. Inst. u. Inst. Angew. Math. Univ. Berlin 4, 
211 (1939). Z 21, 343 

Agnew, R. P., Global versions of the central limit 
theorem, Proce. Nat. Aead. Sei. USA. 40, 800 
(1954). M 16, 268 

Alda, V., On conditional expectations, Czechoslovak 
Math. J. 5, 503 (1955). M 18, 241 

Andersen, E. S., On the number of positive sums of 
random variables, Skand. Aktuarietidskr. 32, 27 
(1949). MI 11, 256 

Andersen, E. 5., On the frequency of positive partial 
sums of a series of random variables, Mat. Tidsskr. 
B. 1950, 33 (1950). M 12, 619 

Andersen, E. S., On the fluctuations of the sums of 


Author: The author’s surname is followed by random variables, Math. Scand. 1, 263 (1953). 
initials only. In the case of complicated surnames, M 15, 444 
| we have used the first capitalized word given in the | Andersen, E.S., On sums of symmetrically depend- 
reviewing journal. Multiple authorship is denoted ent random variables, Skand. <Aktuarietidskr. 





by the symbol @ preceding the surname. The 
journal reference appears after each author’s name 
but the title of the paper is given with the first author 
only. 

Title: This is given exactly as in the reviewing 
journal. Titles of separately bound publications 
(books, reports, ete.) are italicized, and are followed 
by the publisher’s name and address. 

Reference to literature: The name of the journal in 
italies, the number of the volume in bold face, the 


iJ. Research NBS G4B, 55, 


ou (1480 


95237 


0 oo 4 


36, 123 (1953). M 15, 634 
Anis, A. A., On the distribution of the range of partial 
sums of independent random variables, Proce. 
Math. Phys. Soe. Egypt 5, 83 (1953). M 16, 267 
Arfwedson, G., Research in collective risk theory. I, 
Skand. Aktuarietidskr 37, 191 (1954). M 19, 275 


Bachelier, L., Les lois des grands nombres du caleul 
des probabilités (Gauthier-Villars, Paris, 1937). 

Z 16, 170 

Bahn, R.. Uber den Grenzwert der Wahrschein- 

lichkeiten seltener Ereignisse, Deutsche Math. 2, 

698 (1937). Z 18, 33 
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Ballarin, 5., Espressione rigorosa dello searto medi- 
ano nel problema delle prove ripetute nello 
schema di Bernoulli, Vem. Soe. Astr. Ital. 19, 63 


(1948). MI 9, 450 


Banerjee, ID). - Note on the linnit of correlation and 
regression coefficients in mingled records, \Jath. 
Student 9, 155 (1941) M 4, 104 


Barricelli, N. A., L’intégrale relative d’une fonection- 
nelle et ses applications dans la théorie de la dis- 
tribution de probabilité (une courbe, Arch. Math 
Naturvid. 49, 35 (1947). \I 9, 178 

Bavli, G. M., Eine Verallgemeinerung des Poisson- 
schen Grenzwertsitzes, (. PR. Aead. Se7. U RSS 2, 
5OS (1935). Z 12, 216 

Baw ly, G. M., Uber einige Verallgemeinerungen der 
Grenzwertsitze der Wahrscheinlichkeitsrechnung, 


Pec. Math. Moseou, 1, 917 (1936). Z 16, 127 
Bawly, G., Uber den lokalen Grenzwertsatz der 
Wahrsche inlic hkeitsrechnung, Rer. Fae. Sei. Univ. 
Istanbul, 2, 79 (1937). Z 16, 311 
Baxter, G., An analogue of the law of the iterated 
logarithm, Proc. Amer. Math. Soc. 6, 177 (1955). 
MI 16, 112s 
@ Bellman, R.. Recurrence times for the Ehrenfest 
model, Pacifie J. Math. 1,179 (1951). MI 13, 566 
Bellman, R., Limit theorems for non-commutative 
operations. 1. Duke Math. J. 21, 491 (1954). 


M 15, 969 

Bergstrom, H., On the central limit theorem, Skand 
Aktuarietidskr. 27, 139 (1944). M 7, 458 
Bergstrom, H., On the central limit theorem in the 
space Pk >1, Skand. Aktuarietidskr. 28, 106 
(1945). M 7, 459 
Bergstrém, H., On the central limit theorem in the 
case of not equally distributed random variables, 


Skand. Aktuarietidskr. 32, 37 (1949). M 11, 25 
Bergstrém, H., On distribution functions with a 
limiting stable distribution function, .irk. \Jat. 2 


M 15, 237 
S.. Sur les sommes de grandeurs aléatoires 
liées de classes (A,N) et (BN), ©. BR. Doklady 
Acad. Sei. URSS 32, 303 (1941). \I 6, SS 
Bernstein, S., Sur le théoréme limite de la théorie des 
probabilités, Bull. Izvestiya Math. Mech. Inst 
Univ. Tomsk 3, 174 (1946). \I 8S, 471 
Bernstein, S. N. (Editor), The Seient tific Leygae yY of 
P. L. Cebyser. First Part: Mathematics. (Aca- 
demiva Nauk SSSR = Moseow-Leningrad, 1945). 
\l 7. oon 

Retour au probléme de Pévaluation 
formule limite de La- 


$63 (1953). 
Bernstein, 


Bernstein, S. N., 
de Vapproximation de la 


place, Bull. Acad. Sei. URSS, Ser. Math. 7, 3 
(1945). M 5, 41 
Bernstein, S. N., Some remarks concerning the limit 
theorem of Liapunov. Doklady Akad. Nauk 
SSSP 24, 3 (1939). M 1, 340 
Bernstein, S., Nouvelles applications des grandeurs 
aléatoires presqu’indépendantes, Bull. Acad. Se’. 


RSS. Ser. Math. 4, 137 (1940). M 2, 107 
Bernstein, S., Sur une probléme du schéma des urnes 
Q ge variable, ©. R. Doklady Acad. Sei 
RSS 28, 5 (1940). M 2, 229 
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Berry, A. C., The accur: ey of the Gaussian approxi. 
mation to the sum of inde ” ndent variates. Trans. 
Amer. Math. Soe. 49, 122 (1941 \l 2 , 298 
rh irnbaum, Zz. \\., On the properties of il collective 
Amer. J. Math. 62. 787 (1940 \l 2, 106 


Bitvue ‘kov, V.1., A loeal limit theorem for sequences 
of events forming a compound chain of the second 


order, Lerestiya Akad. Nauk SSSR. Ser. Math. 12, 
101 (1948 M 9, 451 
Blackman, J. An extension of the Kolmogoroy 
distribution, dwn. Math. Statist. 27, 513 (1956). 
\I 18, 605 
Blackwell, D.. A renewal theorem, Duke Math. J. 
15, 145 (1948). \l 9, 452 
Blackwell, D., On optimal systems, Ann. Math. 
Statist. 25, 394 (1954 MI 15, 8x2 


@ Blanc-Lapierre, A., 


les fonctions aléatoires statlonnaires 


La loi forte des erands nombres 


pour conti- 
nues, ©. 2. Aead. Sei. Paris 220. 134 (1945). 

M 7, 129 

Blomqvist, N., On an exhaustion process, Shand, 

Ahtuarie tidskr. 35. 20) 1952 \I 14,7 44 71 


almost sure 
1954 
\l 15.7 (2? 


On = convergence of empiric distribu- 


Two theorems on 
Amer. Math. Soe. 


Blum, J. R., 
Proe. 


converg- 


ence, 8, 253 


Blum, J. R., 


tion funetions, ann. Math. Statist. 26, 527 (1955). 
\I 17, 48 
@ Blum, J. R.. A class of stationary processes and 
a central limit theorem, Duke AMMath. J. 24, 
1957). \I 18, 680 
@Blum, J. R.. A class of stationary processes and 
a central limit theorem, Proce. Nat. Aead. Sei. 
S.A. 42, 412 (1956). \I 18, 342 


applicability of the 
Duke Math. J. 12, 43 
Ml 6, 233 

of sums of 
Univ. lé. 

M17, 979 


von Summen 


Conditions of 
numbers, 


Bobroff, A. A., 
strong law of large 
1945). 

Bobrov. A. A., On the 
positive random quantities, 
Jap. ww Mat. 3, 92 (1949 

Bobrov, ,U ber relative Stabilitdt 
oo as ' gufalliger Groében, ( R Acad Se). 
(RSS. 16, 239 (1937 Z 16, 410 

Bobrov, A., A simplified proof of a theorem of A. N. 


relative stabilits 


Voskor. (; IN 


Kolmogorov on the strong law of large numbers, 
( xpeh i Matem. Nauk 2, 194 (1947 \I 9, 519 
Bochner, S., Limit) theorems for homogeneous 
stochastic processes, Proe. Nat. Aead. Sei. USA. 
40, 699 (1954). \I 16, 379 
Bogolvubov, N. N., On certain limiting distributions 


/ che nie 


for sums depending on arbitrary phases, 


Za piski Moskor. Gos. Unir. Fi iha 7 Jae $33 (1945). 
M 7, 314 
Bogolvubov, No N., On the influence of a random 
force on a harmonic oscillator, (chenye Zapiski 
Moskor. Gros. Unie. Fi-ika 77. 5] 1945 ° \l ae 314 
Bohman, H., On a class of orthogonal series, .Ark. 
Mat. 1,13 (1949). M 12, 21 


. 
Borel, Bic. 


Sur les probabilités dénombrables el le 


pari de Paseal, (. R. Aead. Set. «Paris 224, 
1947 \I 8, 280 
Borel, K., Sur les déve ‘loppe ments unitaires normaux, 
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Propagation at oblique incidence over cylindrical obstacles, 
M. P Bachy nski, J. Re search NBS 64D, No 4, 311 (1960). 
Investigations of propagation of short electromagnetic waves 
al oblique incidence over smooth, perfectly conducting 
evlindrical obstacles are described. It is shown that the 
effect of oblique incidence can be considered as a change in 
the effective radius of curvature of the diffracting obstacle. 
The power in the shadow region of a evlindrical 
decreases With angle of obliqueness for horizontally polarized 
waves and can decrease, remain constant, or inerease with 
angle of obliqueness for vertically polarized waves depending 
on the geometry of the propagation link. In all cases, vertical 
polarization gives a stronger field in the shadow region than 
horizontal polarization. In addition it is shown that the 
diffracted field behind an obstruction of uniform radius of 
curvature is the same as that behind an obstacle of uniformly 
varying radius of curvature, provided the effective radius is 
the simne 


obstacle 


Diffraction by smooth conical obstacles, I]. E. J. Neugebauer 
and M. P. Bachynski, J. Research NBS 64D, No. 4,317 (1960) 
Expressions, obtained earlier for the calculation of diffrac- 
tion due to conducting obstacles with smooth cylindrical 
surfaces, are generalized to oblique incidence and to surfaces 
of conical shape. The derivation is based on a generalized 
concept of the Green's function and on the use of corrective 
factors that take the same place as corrections introduced 
by other authors into the theory of diffraction by apertures 
The final expressions for conical obstacles and oblique inci- 
dence are very similar to those for evlindrical obstacles. The 
results are compared with scale model measurements. 


Mode theory and the propagation of ELF radio waves, J. I 
Wait, J. Research NBS GAD, No. 4, 38? (1960). 

The mode theory of propagation of electromagnetic waves 
at extremely-low-frequencies (1.0 to 3000 ©/s) is treated in 
this paper. Starting with the representation of the field as a 


sum of modes, approximate formulas are presented for the 
attenuation and phase constants Certain alternate repre- 
sentations of the individual modes are mentioned. These 


are used as a basis for describing the physical behavior of the 
field at large distances from the source, particularly near the 
antipode of the source At the shorter distances, where the 
range is comparable to the wavelength, the spherical-earth 
mode series is best transformed to a series involving evlin- 
drical wave functions. This latter form is used to evaluate 
the near field behavior of the various field components. 

The effect of the earth's magnetic field is also evaluated using 
a quasi-longitudinal approximation. In general it is indicated 
that if the gvrofrequeney is less than the effective value of 
the collision frequency, the presence of the earth’s magnetic 
field may be negleeted for ELF. When this condition is not 
met the attenuation rrhiay be increased somewhat The 
influence of an inhomogeneous ionosphere is also briefly con- 
sidered and, finally, the propagation of ELF pulses are 
treated It is suggested that certain observed characteristics 
of ELF may be attributed to the 
the current channel in the lightning discharge 


waveforms 


On the diffraction of electromagnetic pulses by curved con- 
ducting surfaces, J. R. Wait and A. M. Conda, Can. J. Phys. 
37, 1384 (1959 

Starting with the known steady-state solutions for diffraction 
by a perfectly conducting convex surface, the corresponding 
transient responses are derived using Fourier-Laplace inver- 
Explicit results are given for an incident wave whict 
Varies With time as a step function 


ston 


inclination of 


Leonard Euler's integral: A historical profile of the gamma 
function, P. J. Davis, Am. Math. Mo. 66, 849 (1959). 

This survey article shows how the gamma function grew in 
concept and in content from the time of Euler to the recent 
treatise of Bourbaki and how in this growth it partook of the 
general development of mathematics over the past two and 
a quarter centuries, 


Confidence intervals for the expectation of a Poisson variable, 
KE. L. Crow and R. 8. Gardner, Biometrika 46, 441 (1959). 

\ table of “optimum” two-sided confidence intervals for the 
mean of a Poisson variable is presented for confidence co- 
efficients 80, 90, 95, 99, and 99.9 percent and all values of 
the variable from 0 through 300. The intervals are compared 
in length with other existing or possible systems of intervals 
for the Poisson mean. The method of calculation is stated, 
and an interesting property of probability sums 
useful in the calculation is derived. 


Poisson 


Use of the equation of hydrostatic equilibrium in determining 
the temperature distribution in the outer solar atmosphere, 
S. R. Pottasch, Astrophys. J.131, No. 1, 68 (1960). 

The temperature distribution from 1.0043 (3000 km) to 20 
solar radii in the sun’s atmosphere is computed from the 
observed density distribution in this region and the assump- 
tion of hydrostatic equilibrium, The temperature distribution 
shows a maximum between 1|.1 and 3 solar radii and a decrease 
in temperature thereafter. This decrease in temperature is 
consistent with Chapman’s suggestion of thermal conduction 
only if loss of energy by radiation is included. Inclusion of a 
radiative energy loss also is shown to invalidate Parker's argu- 
ment against hydrostatic equilibrium out to large distances 
from the sun. 


On the convergence of the Rayleigh quotient iteration for the 
computation of characteristic roots and vectors, VI. (Usual 
Rayleigh quotient for nonlinear elementary divisors), A. M. 
Ostrowski, Arch. Rat. Mech. Anal. 4, No. 2, 1453 (19459). 

In this paper the classical Rayleigh quotient iteration is dis- 
cussed for eigenvalues with non-linear elementary divisors. 
The convergence of the method is only. then satisfactory, if it 
is combined with the accelerating methods of Steffensen and 
Householder, but in this last case it turns out to be at least as 
good as the method of the generalized Rayleigh quotient. 


Tables for the statistical prediction of radio ray bending and 
elevation angle error using surface values of the refractive 
index, B. R. Bean, B. A. Cahoon, and G, D. Thayer, VBS 
Tech Vole Lt (PB151 03 (7960) 5O cents. 


Radio rav bending, vr, and elevation angle error, «, have been 
ealeulated for a wide range of meteorological conditions at 13 
climatically diverse U.S. radiosonde stations. The parameters 
in the observed linear regression equations of 7 and « upon the 
surface value of the refractive index are given for heights of 
0.1 to 70 kilometers and initial elevation angles of the ray from 
0 to 900 milliradians. 


Weighted restricted partitions, M. Newman, Acta Arith. V¥, 


3/1 (1959 


the number of partitions of m into parts not 
divisible by g, and define q,(n by {[Zqi(n)r" n 
In this article recurrence formulas for these coefficients of 
lengths independent of n are derived when q is any of the 


primes, 3. & BD, 4-48 


Let qy(n be 


*=Tq.(n)r 
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A continuous poker game, A. J. Goldman and J. J. Stone, 
Duke Math. J. 27, No. 1, 41 (1960). 

In this paper we derive the solutions of a zero-sum two-person 
poker game in which the players’ hands are independent 
random numbers from the intervals [O, 1]. The game involves 
two bet levels a, b and an ante of 1 unit (a>h>I1 The 
players act alternately, and one of them is permitted a single 
raise, 


Our model subsumes the alternating-bid von Neumann poker 
game of [4] as well as the model solved by Bellman [1]. The 
former arises from the limiting case b=-1, the latter as the 
“equal increments” (Karlin and Restrepo 
[3] have recently solved the equal increments game with n 
rounds of bidding.) The solution exhibits qualitative features 
like those of [4] and [1] and turns out to depend on a decompo- 
sition of [0, 1] into three subintervals corresponding to low 
hands, intermediate hands and high hands. Optimal strategies 
for the players are distinguished among the 
strategies (those which achieve the value of the game against 
every optimal strategy of the opponent) by a specification of 
the average frequeney of bluffing over the range of low hands 
and (for one plaver) by an integral sublinearity condition on 
the frequency of seeing a raise when holding an intermediate 
hand, 


case a—b--hb— 1, 


semioplimal 


List of Titles 


Journal of Research, Section 64A, No. 4, July-August 1960. 
70 cents. 


Gamma _ irradiation of hexafluorobenzene. R. EF. Florin, 
L. A. Wall, and D. W. Brown. 

Behavior of isolated disturbances superimposed on laminar 
flow in a rectangular pipe. Grover C. Sherlin 

Standard of spectral radiance for the region of 0.25 to 2.6 
microns. Ralph Stair, Russell G. Johnston, and BE. W. 
Halbach. 

Photovoltaic effect produced in silicon solar cells by 
gamma ravs. Karl Scharf. 

Phase equilibria in systems involving the rare-earth oxides. 
Part I. Polymorphism of the oxides of the trivalent rare- 
earth ions. R.S. Roth and 8. J. Schneider. 

Phase equilibria in systems involving the rare-earth oxides. 
Part II. Solid state reactions in trivalent rare-earth oxide 
systems. 8. J. Schneider and R. 8. Roth. 

Some observations on the caleium aluminate carbonate hy- 
drates. Elmer T. Carlson and Horace A. Berman 

Acid dissociation constant and related thermodynamic quan- 
tities for triethanolammonium ion in water from O 
50° C Roger G. Bates and Guy F. Allen. 

lonization constants of four dinitrophenols in water at 25° C 
Robert A. Robinson, Marion Maclean Davis, Maya Paabo, 
and Vincent E. Bower. 

Dissociation constant of anisie (p-methoxybenzoic 

the svstem ethanol-water at 25° C 


X- and 


to 


acid in 


Elizabeth E Sager 
and Vincent FE. Bower. 
Preparation of sulfur of high purity. Thomas J. Murphy, 


W. Stanley Clabaugh, and Raleigh Gilchrist 
Tritium-labeled compounds IV. p-Glucose-6-1, b-xvlose-5-f, 
and p-mannitol-/-/. Horace S. Isbell, Harriet L. Frush, 
and Joseph D. Mover. 
Tritium-labeled compounds V. Radioassay of both carbon- 
14 and tritium in films, with a proportional counter 
Horace 8. Isbell, Harriet L. Frush, and Nancy B. Holt. 


Journal of Research, Section 64C, No. 3, July-September 
1960. 75 cents. 


A new method of measuring gage blocks. James B. Saunders 

Gage blocks of superior stability: initial developments in 
materials and measurement. M. R. Meyerson, T. R 
Young, and W. R. Ney. 

Variation of resolving power and type of test pattern 
Kk. Washer and William P. Tayman 


Francis 


A multiple isolated-input network with common output 
C. M. Allred and C. C. Cook. 
Phase angle master standard for 400 cycles per second 


J. H. Park and H. N. Cones. 
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Disturbances due to the motion of a cylinder in a two-layer 
liquid system. Llovd H. Carpenter and Garbis H. Key. 
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Journal of Research, Section 64D, No. 4, July-August 1960, 
70 cents. 


Relation of turbulence theory to ionospheric scatter propa- 
gation experiments. Albert D. Wheelon 

Propagation at oblique incidence over cylindrical obstacles, 
M. P. Bachynski. See above abstracts 

Diffraction by smooth conieal obstacles. H. BE. J. Neuge- 
bauer and M. P. Bachynski See above abstracts.) 

Characteristics of 488 megacveles per second radio signal re. 
flected from the moon. B.C. Blevis and J. H. Chapman, 

The use of polarization fading of satellite signals to study the 
electron content and = irregularities in the ionosphere, 
C. Gordon Little and Robert 8S. Lawrence. 

Note on a test of the equivalence theorem for sporadic BE 
propagation. J. W. Wright and T. N. Gautier 

Davtime attenuation rates in the very low frequeney band 
using atmospherics. W. L. Taylor. 

Measured electrical properties of snow and glacial ice. A. D. 
Watt and E. L. Maxwell 


Half-wave evlindrical antenna in a dissipative medium: ecur- 


rent and impedance. Ronold W. P. King and Charles 
W. Harrison. 

Preface to ELF papers. 

Some ELF phenomena. Ek. T. Pierce 

Mode theory and the propagation of ELF radio waves, 


James R. Wait. See above abstracts.) 


Studies of natural electric and magnetic fields. G. D. Gar- 
land and T. F. Webster 

Natural electromagnetic energy below the ELF range Wal- 
lace H. Campbell. 

Possible application of the system loss concept at ELF, 


Kenneth A. Norton. 

Measurements of the spectrum of radio noise from 50 to 100 
eveles per second. M. Balser and C. A. Wagner. 
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Power loss and operating temperature of tires, R. D. 
Stiehler, M. N. Steel, G. G tichey, J. Mandel, and 
R. H. Hobbs, p. 73 

An indoor tester for measuring tire tread wear, G. G, 
Richey, J. Mandel, and R. D. Stiehler, p. 104 

Measurement of the aging of rubber vulcanizates, J. 


Mandel, F. L. Roth, M. N. Steel, and R. D. Stiehler, 
p. 221 
Standard materials for rubber compounding, F. L 
and R. D. Stiehler, p. 252 

Thermodynamie properties of helium at low temperatures and 
high pressures, D. B. Mann and R. B. Stewart, J. Heat 
Transfer S12, 323 (1959 

Fire research at the National Bureau of Standards, A. F 
Robertson, Fire Research Abstr. Rev. I, No. 4, 159 (1959). 

Refractive indices and transmittances of several optical glasses 
in the infrared, G. W. Cleek, J. J. Villa, and C. H. Hahner, 
tg (pt Am. 49, No. 11, 1090 (1959 

Absolute photometry of the aurora Il. 
emission in the sunlit atmosphere, M. H 
pheric and Terrest Phys 14, 338 (1959). 

Preliminary assessment of the IGY, A. H. Shapley, Proce. Natl. 
Electron. Conf., p. 1 (1958 

Factorial experiments in life testing, M 
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The anomalous inversion in cristobalite, R. F. 
pages (J. 


kinetics of high temperature processes, 228 
Wiley & Sons, New York, N.Y., 1959). 

Digital recording of electrocardiographie data for analysis by 
a digital computer, L. Taback, Ek. Marden, H. L. Mason, 
and H. V. Pipberger, IRE Trans. Med. Electron. ME-6, 
167 (1959). 

Determination of starch in paper. A comparison of the 
TAPPI enzymatic, and spectrophotometric methods, J. L. 
Harvey, B. W. Forshee, and D. G. Fletcher, Tappi 42, 
No. 11, 878 (1959). 

Young’s modulus of various refractory materials as a function 
of temperature, J. B. Wachtman, Jr., and D. G. Lam, Jr., 
J. Am. Ceram. Soc. 42, No. 5, 254 (1959). 

Branched-chain higher sugars. IIT, A 4-C-(hydroxymethvi)- 
pentose, R. Schaffer, J. Am. Chem. Soc. 81, 5452 (1959). 
Vibrational spectrum of cyanate ion in various alkali halide 
lattices, A. Maki and J. C. Decius, J. Chem. Phys. 31, 

No. 3, 772 (1959). 

Glass research at the National Bureau of Standards, C. H. 
Hahner, Glass Ind., p. 588 (1959). 

Applications of the molecular refractivity in radio metrology, 
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resonator, E. Maxwell and F. A. Schmidt, Suppl. 
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1958-1, p. 95 (1958). 
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E. K. Smith, In N. Ath Treaty Organ., Sporadic F 
ionization; Ionospheric Research Meeting, AGARD Avi- 

Panel, Cambridge, England, Sept. 1958 (AGARD- 
ograph 34), p. 129 (1958). 

Low-temperature transport properties of copper and _ its 
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Powell, H. M. Roder, and W. J. Hall, Phys. Rev. 115, 314 
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The double bond isomerization of olefins by hydrogen atoms 
at — 195°, M. D. Scheer and R. Klein, J. Phys. Chem. 63, 
1517 (1959). 
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